Computers and Mathematics with Applications Non-prismatic Timoshenko-like beam model: Numerical solution via isogeometric collocation

shape functions, allows an equal order approximation of all unknown fields, without affecting the stability of the solution. This makes such an approach simple, robust,efficient,andparticularlysuitableforsolvingthesystemofODEsgoverningthenon-prismatic beam problem. Several numerical experiments confirm that the proposed mixed isogeometric collocation method is actually cost-effective and able to attain high accuracy.


Introduction
Non-prismatic structural elements are widely used in several engineering fields, e.g., large span structures, mechanical components, aeronautical applications [1,2]. Unfortunately, their modeling is a non trivial task since the variable geometrical properties of the cross-section lead to scarcely manageable solutions, also when considering simple geometries (i.e., linearly tapered homogeneous beams), loads, and Boundary Conditions (BCs) [3,4]. As a consequence, it is necessary to take advantage of numerical tools (mainly, Finite Element (FE) analysis), which researchers continuously aim at improving in terms of efficiency and effectiveness [5][6][7][8][9]. Furthermore, recent papers have highlighted that smooth changes of the crosssection geometry lead to non trivial stress distributions that deeply influence the whole beam behavior [10][11][12]. In this context, the use of prismatic beam Ordinary Differential Equations (ODEs) simply assuming that the cross-section area and inertia are variable parameters results in a non-effective approach, as noticed since the '60s of the past century [13].
A significant contribution overcoming the problems so far discussed was recently proposed by Balduzzi et al. [14] that derive a novel and effective non-prismatic beam theory based on a simple and consistent definition of kinematics, equilibrium, and constitutive relations. Specifically, all variables are represented with respect to a global Cartesian coordinate system, avoiding possible complications arising from the definition of local coordinate systems. Furthermore, the three cross-section displacements usually adopted in prismatic Timoshenko beam (i.e., cross-section rotation φ, horizontal u and vertical v displacements) and the three internal forces (i.e., bending moment M, horizontal H and vertical V internal forces) are used as unknowns and the resulting beam theory is naturally expressed as a system of six first order ODEs. The numerical results discussed in [14] and the application of the proposed model to several engineering problems [1,4] have already been shown the highly improved performances of the adopted beam theory in terms of ability in predicting stress, strain, displacements, and stiffness.
However, even if it is a 1D theory (hence representing a simplification of the 2D problem), the above mentioned beam model is still represented by a system of ODEs which requires a numerical treatment to obtain solutions in non-trivial cases. Standard numerical approaches such as a displacement-based or mixed FE approach have been already considered in the literature [10], but their use leads to non trivial numerical schemes, for which a high computational cost has to be paid in order to accurately capture the complexity of the solution.
Accordingly, in the present paper, we resort to IsoGeometric Analysis (IGA) (see, e.g., [15] and references therein), and in particular to the family of methods referred to as isogeometric collocation. Initially introduced in [16], aiming at optimizing the computational cost still taking advantage of geometrical flexibility and accuracy typical of IGA, the fundamental idea of isogeometric collocation consists of the direct discretization of the (system of) equations of interest in strong form, adopting the isoparametric paradigm and making use of the higher-continuity properties typical of the spline-based shape functions used in IGA (see, e.g., [17] and references therein).
Isogeometric collocation has been particularly successful in the context of structural elements. In particular, Bernoulli-Euler beam and Kirchhoff plate elements have been proposed by [18], and shear-deformable structural elements have been considered in a number of papers. Mixed formulations both for Timoshenko initially-straight planar beams [19] and for curved spatial rods [20] have been proposed and studied, and then successfully extended to the geometrically nonlinear case [21]. Isogeometric collocation has been moreover successfully applied to the solution of Reissner-Mindlin plate problems in [22], and a new single-parameter formulation for shear-deformable beams, recently introduced by Kiendl et al. [23], has been solved also via IGA collocation.
In the case of the beam theory considered in the present paper, isogeometric collocation seems to be an ideal solution strategy, allowing the creation of a unique combination of a simple and effective analysis tool with a simple and effective structural theory. Therefore, also supported by a number of convincing numerical experiments, we believe that the approach herein proposed may have a significant impact on non-prismatic structural elements modeling.
The paper is structured as follows: In Section 2, we summarize the differential equations governing the beam model under investigation and provide some highlights on the beam model properties. In Section 3, we introduce B-splines and describe the isogeometric collocation method adopted for the numerical solution. In Section 4, we provide several numerical results showing the good behavior of the proposed approach. Finally, in Section 5, we draw our conclusions and discuss possible future developments.

Synopsis of the non-prismatic Timoshenko-like beam
This section recaps the Timoshenko-like beam model ODEs derived in [14]. Readers may refer to the so far mentioned paper for further details on the derivation procedure. The beam behaves under the hypothesis of small displacements and plane stress state and is made of a homogeneous and linear-elastic material. Being l the beam length, we introduce the beam longitudinal axis that coincides with the domain of the functions we are going to introduce. Thereafter, we introduce the beam center line c : L → R and the cross-section height h : L → R + (where R + indicates strictly positive real values). We assume that both the center line and the cross-section height are sufficiently smooth function, as specified in the following. The beam lower and upper boundaries b l , b u : L → R are defined as and, the 2D region occupied by the beam body Ω is defined as ) -3 Fig. 1. Beam model domain, 2D beam geometry, coordinate system, dimensions and adopted notations. Fig. 1 represents the beam longitudinal axis L, the 2D region occupied by the beam body Ω, the Cartesian coordinate system Oxy, the lower and upper boundaries y = b l (x) and y = b u (x), and the center line y = c (x).
Assuming that the lower and upper boundaries b l (x) and b u (x) are unloaded, being σ : Ω → R 2×2 s the 2D symmetric stress tensor and n =  n x , n y  T : b l ∪ b u → R 2 the outward unit vector defined as where the notation (·) ′ indicates the first derivative with respect to the independent variable x and b (x) represents either b l (x) or b u (x), depending on the point where the boundary equilibrium is evaluated. The equilibrium on lower and upper boundaries (σ · n)| b l ∪b u = 0 allows to express the shear stress τ as a function of the axial stress σ x Remark 2.1. Eq. (5) highlights that the shear at the boundary is defined only if the slopes of the lower and upper boundaries are limited, i.e., b ′ l (x) and b ′ u (x) < ∞. Furthermore, discontinuities in b ′ l (x) and b ′ u (x) lead to meaningless discontinuous shear. As a consequence, requiring {c (x) , h (x)} ∈ C 1 (L) is a reasonable hypothesis to ensure that the ODEs are a well-posed problem and the solution has a meaningful physical meaning.
The displacement field s : Ω → R 2 is approximated as in prismatic Timoshenko beam i.e., assuming that the crosssection behaves as a rigid body The beam compatibility is expressed through the following ODEs where the horizontal deformation ε 0 : L → R, the curvature χ : L → R, and the shear deformation γ : L → R represent the generalized deformations.
According to a consolidated approach in structural mechanics [24, Section 1.4], we introduce the internal forces, i.e., the horizontal internal force H : L → R, the resulting bending moment M : L → R, and the vertical internal force V : L → R, respectively defined as follows Being q, m, p : L → R the horizontal, bending, and vertical loads, respectively, the beam equilibrium reads As usual in elementary beam models and taking into account the changes of geometry, we assume that the horizontal stress has a linear distribution within the cross-section and the shear stress has a quadratic distribution (with non-vanishing values at boundaries). In more detail, the distribution of the horizontal stress turns out to be consistent with the cross-section kinematics (see Eq. (6)), whereas the distribution of the shear stress was constructed according to the Jourawsky formula, opportunely modified in order to take into account the boundary equilibrium (5) [25,Section 3.4]. Therefore, the stresses distribution within the 2D region occupied by the body Ω can be expressed as Remark 2.2. The recovery relation (10b) leads to a non-trivial distribution of the shear stress within the cross-section and introduces non-trivial dependencies of the shear stress on horizontal internal force and bending moment [14,25]. This issue is well known since the first half of the 20th century. In particular, referring to the analytical solution of the 2D equilibrium equation for an infinite long wedge, stress distribution can be expressed as the combination of some trigonometric functions [26][27][28]. Furthermore, Timoshenko and Goodier [28] state that a parabolic shear distribution is a reasonably accurate approximation in the case of small boundary slope (i.e., b ′ l (x) and b ′ u (x) should be smaller than tan (15°) ≈ 0.25). Later on, Krahula [29] extends the so far mentioned results to tapered beams, recovering equations substantially identical to (10). Similar results were obtained also by Bleich [30] through an accurate generalization of the Jourawsky theory.
We consider the following simplified expression of stress potential: where E and G denote Young's and shear modulus, respectively. Substituting the stress recovery relations (10) into Eq. (12) and considering the integral with respect to y of the partial derivatives we finally obtain the beam constitutive relations consistent with the so far introduced kinematics and internal forces: G. Balduzzi Remark 2.3. Differently from prismatic beams in which each generalized deformation depends only on the corresponding internal force, constitutive relations (14) highlight that all the generalized deformations depend on all internal forces. In authors' knowledge, the only papers that notice the non-trivial dependence of generalized deformations on all internal forces are [31,10,11,32]. Specifically, deriving a model for a tapered beam, Rubin [32] and Aminbaghai and Binder [31] consider the fact that the bending moment produces also shear deformation and the shear force produces a curvature. Nevertheless both papers do not consider the effect of a non-horizontal centerline and use different coefficients with respect to the ones that naturally result in Eqs. (14). Conversely, Auricchio et al. [10] and Beltempo et al. [11] consider implicitly these aspects through enhanced model derivations.
Substituting the constitutive relations (14) into the compatibility Eqs. (7) we obtain the beam model ODEs in the following mixed form: Obviously, the six first order ODEs (16) must be equipped with six suitable BCs (three per each beam end) in order to ensure the uniqueness of the solution. In Section 4, BCs for every considered numerical example are specified. (1) and highlighted with a bold line in Fig. 1. This modeling strategy is substantially different from the one usually adopted for curved beams that are formulated with respect to a curvilinear coordinate running along the beam centerline [33,34], but has the following main advantages: (i) it avoids complications coming from the handling of a curvilinear coordinate (e.g., curve parameterization, mapping, total derivatives); (ii) it allows to obtain linear and simple ODEs (see Remark 2.5); and (iii) in numerical discretization, it reduces the mapping between real and reference elements to a trivial scaling function. Remark 2.6. As usually done by standard beam theories and according to the adopted kinematics and stress representation, the introduced beam model cannot tackle boundary effects like stress concentrations occurring nearby the constrains [36]. Furthermore, according to Remark 2.2, the formulas for the recovery of stress distributions (10) and consequently the whole beam model provide meaningful solutions only for values of b ′ l (x) and b ′ u (x) smaller than tan (15°) ≈ 0.25.

Numerical solution via isogeometric collocation methods
The aim of this section is to present a short description of B-splines and of their use in isogeometric analysis, followed by an introduction on isogeometric collocation in 1D. We finally present here the application of isogeometric collocation in mixed form for the solution of the system of ODEs describing the beam model introduced in the previous section. The good behavior of the proposed method is then confirmed by the numerical experiments reported in the following section.

Basics on B-splines and their use in isogeometric analysis
We limit the description in the following to the basic concepts on B-splines and isogeometric analysis essential for the proposed method. For more details, the interested reader is referred to [15] and to the references therein.
If internal knots are not repeated, B-spline basis functions are C p−1 -continuous (as in the case of Fig. 2). If a knot has multiplicity k, the basis is C p−k -continuous at that knot. In particular, when a knot has multiplicity p, the basis is C 0 and interpolates the control point at that location. We define the spline space spanned by the basis functions N i,p (ξ ). Following the isogeometric approach [15], the space of B-spline functions defined on a generic interval [a, b] is given by: The images of the knots through the function F naturally define a partition of the interval (a, b), which is referred to as ''mesh'' and is defined by a mesh-size h (being h the largest size of the elements in the mesh).

The adopted mixed isogeometric collocation method
As anticipated in Section 1, isogeometric collocation methods have been introduced in [16] for the numerical solution of ordinary and partial differential equations in strong form by means of an isoparametric approach based on spline functions. They have been shown to be an interesting and efficient alternative to isogeometric Galerkin methods, attracting in particular a lot of attention in the context of structural elements and mixed methods, where equal order approximations seem to be stable and accurate. Such a unique feature can be analytically proven in the univariate context and has been exploited to develop effective Timoshenko-like beam formulations in [20,19].
Following these works, we herein propose to discretize the system of ODEs of Eq. (16) using a mixed isogeometric collocation approach where splines of the same degree are adopted for the discretization of each unknown field. Furthermore, the simple assumptions behind the beam model do not require special treatment of the problem's data (see Remark 2.4). As a consequence, due to the relative simplicity of the numerical examples we are going to discuss, we decide to handle the center line c (x), the cross-section height h (x), and the distributed loads q (x), p (x), and m (x) as analytical functions. Nevertheless, according to Remark 2.1, splines represent the most natural and effective choice for the approximation of both the geometry and the loads in more complex situations. We therefore substitute in Eq. (16) the approximated fields: where ψ(x) indicates one of the six generic unknown functions (i.e., H(x), V (x), M(x), φ(x), v(x), u(x)) andψ j (j = 1, . . . , n) the n corresponding control variables (i.e., the unknown coefficients), while N j,p (x) (j = 1, . . . , n) are the n basis functions of degree p generated by the knot vector Ξ = {ξ 1 = 0, . . . , ξ n+p+1 = 1}. The obtained strong-form equations are then collocated at the imagesx i = F (ξ i ) of the Greville abscissae of the first-derivative spline spaceξ i (i = 1, . . . , n − 1), which can be easily computed as: The resulting algebraic system of equations, consisting of 6(n − 1) equations in the 6n control variables (i.e., n control variables for each of the six unknown functions), has to be finally completed by six suitable BCs to be imposed strongly as additional equations.
Remark 3.1. The possibility of using the same approximation space for all unknown fields, without affecting the stability of the solution, makes the proposed approach particularly simple and effective, and its implementation straightforward. On the other hand, building an equivalent Galerkin approach, dealing with the possible locking effects typical of mixed FEs, may be cumbersome.

Numerical results
The proposed mixed isogeometric collocation method for the solution of non prismatic Timoshenko-like beams is tested on three different examples: (i) a shear-loaded tapered cantilever beam, (ii) an axially-loaded arch-shaped cantilever beam, and (iii) an arch-shaped beam clamped at both the extremities and subjected to self-weight load. Examples (i) and (ii) are taken from [10,14], respectively, which provide also the analytical solution for the ODEs reported in Eq. (16). For Example (iii) we refer to Beltempo et al. [11], where only a numerical estimation of the solution is given due to the high problem complexity, which does not allow to obtain an explicit form of the analytical solution. All investigated examples highlight the non-trivial behavior of a non-prismatic beam body.
For Examples (i) and (ii), the results obtained with the proposed mixed collocation method are reported in terms of convergence analysis. In particular, considering the numerical approximation of a generic solution variable ψ num (x) and its analytical counterpart ψ an (x), we introduce the following error definition where the points x i are uniformly distributed within L and a · n is the number of points used for the error evaluation (the results that follow are obtained setting a equal to 10). It is worth noting that Eq. (24) is an estimation of the error evaluated in L 2 norm. For Example (iii), being the analytical solution not available, the obtained isogeometric collocation results are compared with an overkill solution computed with 2D FE analysis.

Example (i): tapered cantilever
We consider the symmetric tapered beam illustrated in Fig. 3 with the following dimensions: l = 10 mm, h (0) = 1 mm, and h (l) = 0.5 mm. We assume E = 10 5 MPa and G = 4 · 10 4 MPa as material parameters. The beam section at x = 0 is clamped and a vertical load P = −1 N acts on the section at x = l (as illustrated in Fig. 3), x being the axis coordinate. ) -  In the present example, c (x) = 0, therefore, following Eqs. (15) As a consequence, the axial and shear-bending equations become independent. In particular, the ODEs governing the problem read ODEs (25a) have a trivial solution H (x) = u (x) = 0 that the proposed isogeometric collocation method is able to compute exactly (i.e., H num (x) = u num (x) = 0). As a consequence, the error estimations for the horizontal internal force H and for the centerline horizontal displacement u (i.e., err L2 H and err L2 u , respectively) are meaningless and, in the following, we do not consider them.
The analytical solution (also reported in [14]) consists in a constant vertical internal force V (x) = 1 and a linear bending moment M (x) = x − 10. Since both the distributions can be exactly represented by any non-trivial (at least linear) approximating B-spline, with the developed collocation method they can be computed up to machine precision, independently from the number of used collocation points n. As a consequence, their convergence is trivial and it is not reported in the following. Fig. 4 reports the variation of the error estimations err L2 v and err L2 φ for different number of collocation points n, chosen always equal for each field. All polynomial orders are chosen equal to p = 2, 4, and 6, respectively, for the shape functions of each of the six independent fields. The obtained results show that the convergence rates of the collocation method agree with the expected ones [19].

Example (ii): Arch-shaped cantilever beam
We consider the arch-shaped beam depicted in Fig. 5 with the following dimensions: l = 10 mm, ∆ = 0.5 mm, h (l) = h (0) = 0.6 mm, q = 1 N/mm. We assume E = 10 5 MPa and G = 4 · 10 4 MPa as material parameters. For this example, the beam centerline and the section height are defined as   respectively. The beam section at x = 0 is clamped and a horizontal load Q = 0.6 N acts on the centroid of the section at x = l (as illustrated in Fig. 5). Consequently, the BCs for the ODEs (16) read The analytical solution is provided in [14] and consists in a constant horizontal internal force H (x) = 0.6, in a quadratic distribution of the bending moment M (x) = 0.6 · c (x), while the vertical internal force vanishes (i.e., V (x) = 0).
Since, also in this case, the distributions of the internal forces can be exactly represented by any (at least quadratic) approximating B-spline, with the proposed isogeometric collocation method they can be computed up to machine precision, independently from the number of adopted collocation points n. Consequently, in the following we do not report convergence results of internal forces solution.
We rather represent in Fig. 6 the error estimations err L2 v , err L2 u , and err L2 φ for each cross-section displacement evaluated using different number of collocation points n. Also in this case, we keep equal the number of collocation points for each field, as well as the polynomial orders (i.e., p = 2, 4, and 6, respectively) of the shape functions of the six fields. Expected convergence rates, in agreement with Beirao da Veiga et al. [19], are obtained.
As shown in Fig. 7, the beam is clamped at both extreme cross sections and subjected to self-weight load p equal to the artificial value of 50 kN/m 2 . As a consequence, the vertical load is defined as p (x) = −50 h (x) and the BCs read In order to provide a reference solution we model the 2D beam body with 2D FEs, using the commercial software Abaqus (Dassault Systémes, Paris, France). Specifically, we use quadrilateral bilinear elements with a characteristic length of 0.001 m leading to a structured (overkill) mesh of 6 513 396 elements. Fig. 8 reports the numerical solution for the problem under investigation using the proposed IGA collocation method with n = 8 and p = 4 for each field. The reference FE solution is reported at 21 sections, equally distributed along the beam length.
It is worth noting that, already for n = 8, an accurate numerical solution is attained in terms of both internal forces and displacements, with negligible computational costs as compared with 2D FE analysis. Finally, Eqs. (10) allow to accurately recover the cross-section distribution of stresses (as extensively discussed in [11,14]) and therefore to perform accurate analysis of the cross-section strength. Therefore, Fig. 8 demonstrates that the herein proposed method is a promising analysis tool allowing to estimate all the quantities of interest for practitioners also for complex structural elements with a computational effort similar to that required by standard 1D FEs.

Conclusions
This paper presents an innovative numerical approach for the solution of the ODEs governing the problem of planar nonprismatic Timoshenko-like beams that naturally gives rise to a mixed formulation. The proposed scheme is based on mixed isogeometric collocation, where B-splines of the same order are used for the approximation of all (six) unknown variables, and the governing equations of the problem are solved in strong form. The method is simple, cost-effective, and robust. Finally, several numerical experiments have confirmed its good accuracy, even in the case of coarse meshes. If compared with FEs, the most diffused numerical tools for structural analysis, the herein considered numerical approach presents the following main advantages: (i) it avoids the possible locking effects typical of mixed FEs resulting to be more robust and simple to implement, and (ii) it requires a computational effort smaller than (at most similar to) standard 1D FEs.
Further developments will include the generalization of the proposed combination of the beam model and the collocation method to more complex situations like multilayered and 3D beams.