Large nonuniform bending of beams with compressible stored energy functions of polynomial-type

The large bending of beams made with complex materials nds application in many emerging elds. To describe the nonlinear behavior of these complex materials such as rubbers, polymers and biological tissues, stored energy functions of polynomial-type are commonly used. Using polyconvex and compressible stored energy functions of polynomial-type, in the present paper the equilibrium problem of slender beams in the fully nonlinear context of nite elasticity is formulated. In the analysis, the bending is considered nonuniform, the complete three-dimensional kinematics of the beam is taken into account and both deformation and displacement elds are deemed large. The governing equations take the form of a coupled system of three equations in integral form, which is solved numerically through an iterative procedure. Explicit formulae for displacements, stretches and stresses in every point of the beam, following both Lagrangian and Eulerian descriptions, are derived. By way of example, a complete analysis has been performed for the Euler beam.


Introduction
Today we attend to an extraordinary evolution of technical applications based on the nite bending of beams. For example, electric sensors have recently been applied to human motion monitoring and human-machine interaction [1], [2]. These sensors are integrated into exible beams, oering some advantages such as extreme durability and high reproducibility [3]. Piezoelectric actuators are used for their excellent guiding accuracy during bending of beams. These devices are capable of converting small change in length into a large vertical displacement [4]. Large bending actuators are also made using shape memory alloys (SMA). In this case, a contractile wire, usually in nickel-titanium, is applied on the surface of a exible beam or strip [5]. Robots based on extremely compliant components constitute an emerging eld (so called soft robots) [6].
These high-tech devices can be stimulated to move by changes, for example, in the intensity of the light [12], environmental humidity [13] or electric eld [14], etc.
The analysis of bending beams in the context of nite elasticity is considerably complicated by the intrinsic nonlinearities in the mathematical formulation due to the kinematics, the constitutive law and the equilibrium conditions. The main contributions to formulation and resolution of this equilibrium problem are briey outlined in the following.
The large uniform bending of a plate was rst studied by Seth [15]. Based on the semi-inverse approach, he assumed the deformed conguration of the plate like a short circular cylindrical shell, keeping valid the Bernoulli-Navier hypothesis for cross sections. In addition, he assumed a plane strain condition along the axis of the cylinder and a linear constitutive relationship.
The uniform exion of an elastic block has been investigated by Rivlin [16].
The bending transforms the block into a cylinder with the base in the shape of a circular crown sector. No displacements along the axis of the cylinder have been considered, making the problem as a matter of fact two-dimensional. Rivlin also showed that, in the case of a neo-Hookean material, the surface tractions necessary to induce the hypothesized deformed conguration are statically equivalent to two equal and opposite couples acting at the end faces. Several papers based on Rivlin's solution can be found in Literature (see, e.g., [17], [18], [19] and [20]).
A closed-form solution of a compressible block made of Hencky material has been obtained by Bruhns et al. [21], giving explicit relationships for the bending angle and bending moment as functions of the circumferential stretch. Haughton [22] considered Neo-Hookean and Ogden materials, while Triantafyllidis [23] assumed incrementally linear constitutive relationships to reproduce un elasticplastic behavior. A similar study was carried out by Coman and Destrade [24] and then extended to layered solids by Roccabianca et al. [25].
By passing from a block to a beam a wide amount of studies dealing with large deections can be found in the Literature for this special structural element under several loading and clamping conditions. A signicant part of these studies is based on the solution of the Elastica according to the well-known Euler-Bernoulli law for bending (cf. Love [26]). Shield [27] investigated the problem of a beam under pure bending by assuming large displacements but small strains. In particular, he retrieves the Lamb's solution [28] for the deec-tion of the middle surface of the beam.
After some numerical studies about a cantilever beam subjected at its free edge to a concentrated vertical load, Wang et al. [29] proposed a straightforward numerical method based on the nite dierences to solve the equilibrium problem of beams under dierent load distributions. A comprehensive review on the application of the Elastica can be found in the book by Frisch-Fay [30].
Later, the same subject has been reconsidered by Wang [31], conrming the validity of the Newton-Raphson numerical method in solving the transcendental equations governing the problem, and by Holden [32], who solved the problem through a fourth-order Runge-Kutta procedure. However, in these works, a linear relationship between the curvature and the bending moment has been adopted.
In the framework of Finite Element Methods (FEM) for nonlinear analysis of structures, many works concerning the large displacements and large rotations of beams have been carried out. As an example, Bathe and Bolourchi [33] reported a Lagrangian formulation based on incremental equilibrium equations and a proper decomposition of stresses and strains. Cubic interpolating functions are assumed to describe the displacement eld related to bending. A straightforward parametrization of the equation of motion suitable for FEM has been proposed by Simo [34]. In [34], the deformed conguration of a beam is completely described by an orthogonal matrix, from which both the rigid rotations of cross sections and the position of the centroids can be inferred. In addition, in this work, it is shown that the formulation reported by Reissner [35] is exactly retrieved when a plane problem is considered. After the paper by Reissner [36], a lot of studies about the problem of beams under nite displacements have been carried out (e.g., [37], [38], [39], [40], [41] and [42]). However, in these works, neglecting the quadratic part of the Green-de Saint Venant strain tensor, small strains are considered and a linear constitutive relationship is adopted.
In all the above quoted works relating to the bending of blocks and beams under nite displacements, the pure deformation of cross sections has been completely neglected. In this way, the modeling of the problem is simplied substantially, since the displacement eld is assumed to be plane, renouncing to describe a phenomenon which in reality is purely three-dimensional. Specically, the transverse deformation, always coupled with the longitudinal inexion of a solid, is known as anticlastic eect.
Recently a realistic structural model has been proposed for the large uniform bending of beams by Lanzoni and Tarantino [43]. Following a semi-inverse approach, a three-dimensional kinematic model, where the longitudinal bending is accompanied by the transversal anticlastic deformation of cross sections, has been formulated. The theoretical model proposed in [43] has duly been veried through both numerical analyses and experimental campaigns. The numerical analyses are based on FEM simulations, whereas a test equipment prototype has been designed and manufactured for the large bending of slender beams [44], [45], [46] and [47]. In general, numerical and experimental analyses have provided results that are consistent with those given by the theoretical analysis [48].
In view of the many technical applications that can be conducted with this new structural model, the need to generalize the adopted constitutive law has emerged. With this regard, in the present paper the equilibrium problem for beams in the context of nite elasticity has been completely reformulated for a wide class of materials whose stored energy function takes a polynomial expression. For this class of materials, the stored energy function is expressed as a power sum of the principal stretches which satises the polyconvexity conditions. In Literature, this density function has been extensively applied to model the nonlinear behavior of complex materials such as rubberlike solids, polymers and biological tissues [49], [50] and [51].
The paper is organized as follows. The kinematic model of inexed beams with variable curvature is recalled in Section 2. Using a semi-inverse approach, the displacement eld derived from the kinematic model contains four unknown functions, which are evaluated in Section 3 by imposing equilibrium conditions.
The governing equations constitute a coupled system of nonlinear equations in integral form, which is solved numerically through an iterative procedure.
Considering the Euler beam, a numerical application has been performed in of the three-dimensional Euclidean space E. The symbols B, H and L respectively denote the width, height and length of the body. As is typical in the case of beams, the length L is predominant on the both transverse dimensions B and H. As shown in Fig. 1, three-dimensional beams vertically inexed are considered. Although the formulation will be developed for beams with a rectangular cross section, it can readily extended to beams with a generic cross section provided that the symmetry with respect to Y axis is maintained.
The undeformed congurationB of the beam is assumed as the reference conguration, whereas the deformed conguration is given by the deformation 1 It should be kept in mind that the internal constraint of incompressibility, especially in the case of large deformations, can signicantly aect the shape assumed by an inexed body [52]. f :B → V, 2 that is a smooth enough, injective and orientation-preserving (in the sense that det Df > 0) vector eld. The deformation of a generic material point P belonging toB can be expressed by the well-known relationship f (P ) = s(P ) + id(P ), where id(P ) and s(P ) = u(P )i + v (P )j + w (P )k (2) are the position and displacement vectors of the point P. In the vectorial equation (2), functions u( P), v(P) and w(P) are the scalar components of s(P), whereas i, j and k are the unit vectors. In the sequel, an Eulerian coordinate system {O, x, y, z} is also used (cf. Fig. 2).
To derive the displacement eld of an inexed beam under variable curvature a semi-inverse approach, which involves the denition of a kinematic model, is followed. This model is based on the following three basic hypotheses.
1. In the bending of the beam, cross sections maintain their planarity (Bernoulli-Navier hypothesis). That is, plane cross sections, orthogonal to Z axis (cf. Fig. 1), remain as such after the beam has been inexed.
2. Due to the longitudinal inexion also cross sections are transversely inexed, but with opposite curvature (anticlastic eect). This transversal inexion is assumed to occur with constant curvature (cf. Fig. 2c).
3. Slender beams with compact cross sections are considered.
2 V is the vector space associated with the Euclidean space E. The rst assumption, which is very popular in the linear mechanics of beams under pure bending, predicts the conservation of the planarity of cross sections after the deformation. For beams that satisfy the third assumption, the reliability of such hypothesis in nonlinear theory has recently been checked both from the experimental and numerical points of view in [44], [45] and [48]. The second hypothesis also becomes acceptable when beams satisfy the third hypothesis [48].
Based on the three previous hypotheses, the displacement eld is derived as the outcome of the coupled eects generated by both longitudinal and transversal curvatures [43]: where v N (s) = v(0, 0, Z ) and w N (s) = w(0, 0, Z ) are the displacement components of the centroid of the generic cross section, whereas θ N (s) = θ(0, 0, Z ) represents its rotation. The function r(s) denotes the anticlastic radius (cf. Fig.   2c). As shown in [48], in the case of slender beams, the beam axis does not undergo variations in length (λ Z = 1). Thus, although the beam is not, its axis can be considered inextensible. Consequently, Z = s, where s is the curvilinear abscissa measured in the deformed conguration (cf. Fig. 2a).
Using (1) and (3), the components of the deformation gradient F can be expressed as 3 Given the polar decomposition theorem, F = RU, it is immediate to write the deformation gradient (4) as product of the rotation tensor R by the stretch tensor U, where 4 On the basis of the hypotheses formulated for the kinematic model, the deformation gradient (4) describes the deformation state of each single point of the beam. The rotation tensor (5) highlights that each point undergoes a longitudinal rotation θ N (s), which is the same for the entire cross section, and a transversal rotation β(s) = X/r(s), within the cross section, that depends on the variable X (cf. Fig. 2c). The stretch tensor (6) characterizes the pure deformations. It is diagonal because the reference system {O, X, Y, Z} is principal for the deformation state derived from (3). The principal stretches λ X , λ Y and λ Z in (6) have the following expressions: being In the previous formula, R(s) denotes the longitudinal radius of curvature of the beam axis (cf. Fig. 2b).
On the basis of the hypotheses formulated at the beginning of this Section, the nonlinear displacement and deformation elds for beams subject to variable bending moment have been recalled. In these elds there are four functions that have not yet been determined: The three functions related to the deformed axis of the beam, v N (s), w N (s) and θ N (s), and the anticlastic radius, r(s). These unknown functions will be evaluated in the next Section by adopting a specic constitutive law and subsequently by imposing the equilibrium conditions.

Constitutive relationships and equilibrium
Constitutive properties of a hyperelastic material are described by the stored energy function ω. If the function ω is frame-indierent, homogeneous and isotropic, then it depends only on the principal invariants I i , with i = 1, 2 and 3, of the left Cauchy-Green strain tensor B = FF T (which coincide with those of the right Cauchy-Green strain tensor C = F T F) 5 With these assumptions, the constitutive law (T R = ∂ω/∂F) takes the following form: where the tensor T R denotes the (rst) Piola-Kirchho stress tensor. Being BF = RU 3 and F −T = RU −1 , the tensorial equation (9) can be rewritten as 5 The following notations: A = trA T A 1/2 for the tensor norm in the linear tensor space Lin and A = (detA)A −T for the cofactor of the tensor A (if A is invertible) are used.
From the kinematic model, the tensors R and U are known (see (5) and (6)).
Being the tensor U diagonal also the tensor S in (10) is such Equilibrium requires that the following vectorial equation be satised locally: In the absence of body forces b and computing the scalar components of the material divergence of T R , the vectorial equation (11) provides the following system composed of three partial dierential equations: The derivatives S J,J = ∂S J ∂J , for J = X, Y, Z (no sum) are reported in the Appendix. System (12), governing equilibrium conditions locally, shows an exceptional complexity and, consequently, a practical impossibility of being solved formally. But above all, it must be taken in mind that, since the shape of the displacement eld has been assigned a priori through the hypothesized kinematic model, system (12) cannot provide the exact equilibrium solution for all internal points of the beam.
Nevertheless, some correct information can be obtained by imposing the equilibrium in special points of the beam. In particular, the points of the beam axis show a kinematics completely and properly described by the two radii R(s) and r(s). 6 With the aim of rewriting the system (12) for the points of the beam axis it is observed that for them 6 The radius r(s) in the kinematical model is assessed at the points belonging to the beam axis and then it has been employed to describe the displacement of all other points of the cross section. Therefore, as one moves away from the beam axis, the value of r(s) can become approximate. 7 The three S J terms are the same because, at the beam axis, stretches and invariants are the same.
The expressions of S X,X and S Z,Z are similar because, on the basis of the kinematic model, it results: λ X,X = I 1,X = I 3,X = 0 and λ Z,Z = I 1,Z = I 3,Z = 0. The expression of S Y,Y is dierent for the dependence on the variable Y of λ Y , I 1 and I 3 .
, λ Z,Z = 0, and system (12) specializes into (β(s) = 0) Requesting that the stress be zero in the absence of deformation (λ X = λ Y = λ Z = 1), it is obtained: S X = S Y = S Z = 0, and system (14) reduces to Since this system must be solved for any θ N (s) values, it can be further reduced System (15) expresses the equilibrium conditions for the points belonging to the beam axis. In particular, the rst equation governs the equilibrium along the X axis. The second and third one, in coupled way, govern the equilibrium along the Y and Z axes.
In order to solve the system (15) some (only sucient) conditions can be immediately identied to have S X,X = S Z,Z = 0 for X = Y = 0 ( 8 ) 8 Obviously these derivatives can be dierent from zero for the other points of the cross section.
The second equation, S Y,Y = 0, provides then an expression that links the anticlastic radius r(s) to the local value of the longitudinal radius of curvature R(s). In general, system (15), as set up, is able to provide only information on the radius r(s). By inserting the expression for r(s) in the kinematic model, the radius R(s) remains to be determined. In the following other equilibrium conditions will be exploited to evaluate the others unknown functions.
Lagrangian stresses were expressed by the Piola-Kirchho stress tensor T R , while the stress measure coherently employed in the spatial conguration is that of Cauchy. The Cauchy stress tensor T is obtained from the Piola-Kirchho stress tensor T R through the well-known transformation which given (4), (5), (6) and (10) where Since Cauchy stresses are Eulerian functions, in the above formulae it is more appropriate to use the Eulerian stretches derived from (7) [43], [48] The principal directions of stress are the eigenvectors associated with these eigenvalues. The principal direction corresponding to the eigenvalue T 3 is the unit vector orthogonal to the plane Ω (cf. Fig. 2c) with components (0, − sin θ N (s), cos θ N (s)). The others two eigenvectors are any two unit vectors orthogonal to each other and belonging to the plane Ω .
To move forward in the problem formulation it is now necessary to dene the constitutive relationship. The use of stored energy function of polynomial-type is probably the most general and popular approach to describe the nonlinear constitutive behavior of complex materials such as rubbers, polymers and biological tissues. In this case, as suggested by Ogden, the rst two terms of the stored energy function can be expressed as sum of powers of principal stretches [49], [50], [51] 9 = δ, is a convex function that satises the growth conditions both as δ → 0 + and δ → ∞. If the following inequalities are fullled: α 1 ≥ 2, 1 α1 + 1 β1 ≤ 1, then, no loss of ellipticity occurs in the equilibrium equations and the existence of solutions of the boundary-value problem is assured in a properly chosen function space [51].
Hereinafter the expression proposed by Ciarlet and Geymonat [53] is chosen for where c and d are positive constants.
Note that from (21) and similarly the stored energy function for compressible neo-Hookean materials A classic constitutive law for hyperelastic materials is usually expressed in terms of derivatives of the stored energy function with respect to the deformation invariants, ω i = ∂ω ∂Ii , for i = 1, 2 and 3 (cf. (9)), while the variables used in the stored energy function (21) are the three principal stretches λ i , for i = 1, 2 and 3. To use (21) it is therefore necessary to establish correlation formulae between the derivatives of ω with respect to the invariants I i , ω i = ∂ω ∂Ii , and 9 The principal stretches λ X , λ Y and λ Z are renamed respectively with the symbols λ 1 , λ 2 and λ 3 . 10 The compressible Mooney-Rivlin stored energy function was used, for example, in [54], [55], [56] and [57]. those with respect to the stretches λ i , ∂ω ∂λi . Using the chain rule, ∂ω ∂λi = ∂ω ∂Ij ∂Ij ∂λi , the following expressions can be written between the two groups of derivatives: By inverting this system the correlation formulae are obtained From (21) the derivatives to be introduced into system (26) are In the kinematics of the beam illustrated in Section 2, the transversal stretches λ 1 and λ 2 are equal. Then, performing the passage to limit λ 1 → λ 2 , system (26), with (27), gives (λ 1 = λ 2 = λ) In the even more particular case in which all three stretches are equal, with a further passage to limit, formulae (28) are transformed into (λ 1 = λ 2 = λ 3 = λ) By introducing (28) in (10), Piola-Kirchho stresses can be computed for each point of the beam using the stored energy function (21).
Among the constitutive constants involved in (21) a relationship can be established by imposing that, in the absence of deformation, the stress vanishes.
By setting β(s) = θ N (s) = 0 into (5), the stresses T R,ij , with i = j, for i, j = 1, 2, 3, are zero, whereas the diagonal component for λ i = 1 are 11 By inserting in this condition the derivatives (29), written for λ i = 1, the following relationship for the evaluation of the constant d is obtained: Using this expression, at the beam axis, it occurs that S X = S Y = S Z = 0.
The equilibrium of the points belonging to the beam axis is governed by system (15) and the derivatives S J,J , with J = X, Y and Z, are specied in the Appendix. In the expressions S J,J there are the partial derivatives of ω i with respect to the variables X, Y and Z, which are evaluated at the points of the beam axis: X = Y = 0 and Z = Z. Therefore, using (7) and (28) and computing partial derivatives, the following nine expressions are obtained: Since only the derivatives with respect to variable Y are dierent from zero, the rst and third equations of system (15) are identically satised S X,X = 2ω 1,X + 4ω 2,X + 2ω 3,X = 0, S Z,Z = 2ω 1,Z + 4ω 2,Z + 2ω 3,Z = 0. 11 A similar position has been used, for example, in and [58], [59] and [60].
As already mentioned, the second equation of system (15) can be used to derive the relationship between the two radii r(s) and R(s). Thus, using (13), (32) and (29) where the constant d can be calculated with (31). If the radius r(s) takes this expression, the system (15) is fullled, although the radius R(s) is still to be determined.
In the special case of a Mooney-Rivlin material, (28), (31) and (33) are reduced respectively to These expressions coincide with those obtained in [43] and [48].
Once the Piola-Kirchho stresses have been evaluated using (21) wherex(s) = ρ sin β,ỹ(s) = r(s) − ρ cos β, and This relationship depends on the unknown radius R(s) and for this reason we will refer to it as the moment-curvature relationship in the general form. In fact it, derived in the fully nonlinear context of nite elasticity, relates the external bending moment m ext x (s) with the longitudinal curvature R(s) −1 , both evaluated along the deformed beam axis. The moment-curvature relationship (36) will be used below to determine the displacement eld of the beam axis.
The two displacement components v N (s) and w N (s) allow evaluating the position of the deformed axis of beam, and they can be written in terms of At this point, the solving procedure for the nonlinear equilibrium of inexed beams can be organized into three successive steps.
The rst step of procedure concerns the resolution of system composed of equations (36) and (37), whose unknowns are the three kinematic functions v N (s), w N (s) and θ N (s), that is the displacements and rotations of the beam axis. This system is coupled and nonlinear. It can be solved through the iterative numerical procedure illustrated in Appendix A of [43]. In each single interaction the longitudinal curvature R(s) −1 is updated using (8) and consequently also the transversal radius r(s) is modied according to (33). Once the iterative procedure has reached convergence, in addition to the three kinematic functions, the bending moment m ext x (s) and the two radii R(s) and r(s) are known for the specic case treated.
The second step deals with the determination of the three-dimensional shape assumed by the beam in the deformed conguration. Since functions v N (s), w N (s), θ N (s) and r(s) are available, the displacements of each point of the beam can be evaluated with (3).
The third step is dedicated to the analysis of stretches and stresses at the points of each cross section of the beam. For this purpose, the Lagrangian expressions (7) and (10) or the Eulerian expressions (20) and (17) can be used.
Following this procedure, the complete solution is obtained. However, it is necessary to keep in mind that the above solution was not obtained through a classic boundary-value problem. But a semi-inverse method has been applied, which hypothesizes the form of the solution by assigning the displacement eld unless four unknown functions. Operating in this way inevitably some approximations will be present with respect to the exact solution, the determination of which, however, remains excessively complex.
To estimate the accuracy of the solution proposed, the vectorial equilibrium equation (11) can be used, checking if it is satised locally, that is for each single point of the beam. For this purpose, the three scalar equations derived from (11) are rewritten in dimensionless form, so that their comparison with the scalar zero takes full meaning. Therefore, by evaluating how much these three equations deviate from zero, it is possible to measure the accuracy of the solution obtained. In [52], it has been shown that, for the points belonging to the beam axis, the equilibrium equations are exactly satised (see also [46] and [48]).
In fact, for these points, the displacement eld is correct being it kinematically

Numerical application
As an example of application, the Euler beam is considered. This classic bifurcation problem will be reformulated, and analyzed in all its aspects, on the basis of the beam model developed in this paper. A beam with the following geometrical dimensionless parameters: H = 1, B = 2, L = 15 is taken into account. The constitutive behavior of this beam is described by (21), where the series are truncated at the third term (L = 2, M = 1) and the following exponents: α 1 = 4, α 2 = 2 and β 1 = 2 and dimensionless constitutive constants: A 1 = 0.1, A 2 = 1, B 1 = 0.1 and c = 1 have been chosen [61]. Adopting these parameters, (31) gives d = 4.8.
As it is well known, immediately after the bifurcation, which occurs when the axial compressive force reaches the Euler load P CR = π 2 EJ X L 2 , the beam assumes one of the two adjacent stable equilibrium congurations. 12 In these symmetrical congurations, the shape of the beam axis is described by the following where v 0 is the vertical displacement at Z = L/2. In the classical analysis of the Euler beam, the parameter v 0 remains undetermined. Otherwise, in the present case study, the three-dimensional shape taken by the beam after bifurcation is completely determined. With reference to the upper deection curves, Fig. 3a shows the equilibrium post-buckling paths of the Euler beam as the multiplier µ of P CR increases. In this gure, the beam 12 The symbol E denotes the Young's modulus and J X the moment of inertia of the undeformed cross section with respect to the X axis. in the dierent deformed congurations, as well as the thicknesses, are plotted exactly in scale. As µ grows, the right sliding support moves towards the left hinge. For µ = 2.184, the support has traveled the full length of the beam, reaching the hinge. In this particular situation, the beam axis takes the shape of a drop. After this value of µ, the support exceeds the hinge, and the beam forms a loop.
To apply the numerical procedure exposed at the end of Section 3, the Euler beam has been discretized into 100 sampling nodes. The rst-step trial solution has been taken proportional to the initial sinusoidal deection curve. Up to µ = 1.5, the convergence towards the equilibrium solution requires less than 15 iterations. For higher values of µ, the number of subdomains has been increased to improve the convergence. In particular, the equilibrium solutions for µ = 2.184 and µ = 3. The longitudinal Cauchy stress T 3 is shown in the Fig. 3c for each point of the cross section by means of contour lines. The stress component T 3 assumes the maximum tensile value at the upper edge, T 3 = 2.867, and that of compression at the lower edge T 3 = −3.510. Note that the neutral line for the stress (T 3 = 0) does not coincide with the neutral axis for the deformation (λ z = 1).
The constitutive constants used in the application have been made dimensionless, dividing them by the constant A 2 , to which the unit value was assigned.
If the elastic constants are without dimensions, then, stresses are also such. Likewise, the geometrical dimensions of the beam have been rendered dimensionless, dividing them by the height H, which was assumed to have a unit value. In the same way, the variables X, Y and Z are normalized. Since these variables and stresses have been normalized, equilibrium equations (11) are also dimensionless and can therefore be compared with scalar zero. Evaluating therefore how much the three scalar equilibrium equations deviate from zero, it is possible to measure the accuracy of the solution obtained for the Euler beam.
Equilibrium equations are exactly satised at points belonging to the beam axis [52]. This is because for these points the displacements eld is correct, being it kinematically compatible and equilibrated (system (15) has indeed been satised). On the other hand, when moving away from the beam axis, the equilibrium equations are no longer so well satised and some approximations appear as the edges of the cross sections are approached, since for these points the hypotheses underlying the kinematic model become less accurate.
In the light of the above, a numerical analysis was carried out, evaluating the equilibrium equations (11) at all points of the most stressed cross section (s = L/2 and µ = 3.3). Results are delivered in the diagrams illustrated by Figs. 3d, 3e and 3f for the equilibrium along the X, Y and Z axis, respectively.
In these diagrams some contour lines, which join the points where the equations (11) give the same numeral values, are shown. As can be seen from Figs. 3d, 3e and 3f, even in the most critical case of the cross section subject to the maximum bending moment, the numerical values are very close to zero at each point of the cross section, showing that the solution obtained is actually accurate.

Conclusions
In the framework of fully nonlinear elasticity, the equilibrium problem of nonuniform bending of beams with compressible stored energy functions of polynomialtype has been formulated. By means of a kinematic model, which also takes into account the anticlastic inexion of the cross sections, the three-dimensional displacement eld of the beam has been described using four unknown func- By assembling the derived formulae, the governing equations are obtained.
These take the form of a coupled system of three equations in integral form, which is solved numerically through an iterative procedure. Solving this system, the displacement eld has been obtained and then the shape assumed by the beam in the deformed conguration, the stretches and stresses in every point of the beam (following both Lagrangian and Eulerian descriptions) have been assessed.
As numerical application, the nite bending of the Euler beam has been considered. Unlike the classic solution, the shape assumed by the deformed beam has been precisely determined and shown in a series of graphs, even in the cases with very high load multipliers. The vertical proles of longitudinal stretch and stress for the deformed cross section, subject to the maximum bending moment, are displayed. Finally, the accuracy of the solution obtained has been estimated by checking a posteriori that the dimensionless equilibrium equations practically do not diverge locally from zero.
Appendix. Expressions of the derivatives S J,J In system (12) the derivatives S J,J = ∂S J ∂J , for J = X, Y, Z (no sum), assume the following forms: