Beam-like models for the analyses of curved, twisted and tapered HAWT blades in large displacements

The continuous effort to better predict the mechanical behavior of wind turbine blades is related to lowering the cost of energy. But new design strategies and the continuous increase in the size and flexibility of modern blades make their aero-elastic modeling ever more challenging. For the structural part, the best compromise between computational efficiency and accuracy can be obtained by schematizing the blades as suitable beam-like elements. This paper addresses the modeling 10 of the mechanical behavior of complex beam-like structures, which are curved, twisted and tapered in their reference state, undergo large displacements, 3D cross-sectional warping and small strains. A suitable model for the problem at hand is proposed. It can be used to analyze large deflections under prescribed loads and allows the 3D strain and stress fields in the structure to be determined. Analytical results obtained by applying the proposed modeling approach are illustrated.


Introduction 15
In the process of improving horizontal axis wind turbines (HAWT) performance new methods are continuously being sought for capturing more energy, developing more reliable structures, and lowering the cost of energy (Wiser, 2016). Such goals can be achieved through the use of advanced materials, the optimization of the aerodynamic and structural behavior of the blades, and the exploitation of load control techniques (see, for example, Ashwill 2010, Bottasso 2012, Stäblein 2017. But new design strategies and the increasing flexibility of those structures make the modeling of their aero-elastic behavior ever 20 more challenging. For the structural part of this modeling, schematizing the baldes through suitable beam-like elements can be the best compromise between computational efficiency and accuracy. But modern blades are complex beam-like structures. They can be curved, twisted and also tapered in their unstressed state. Even not considering the complexities related to the materials and loading conditions, their shape alone is sufficient to make the mathematical description of their mechanical behavior a very challenging task. This paper addresses the modeling of the mechanical behaviour of structures of 25 this kind, with a particular focus on their main geometrical characteristics, such as the twist and taper of their cross-sections, the in-and out-of-plane warping of their cross-sections, and the large deflections of their reference center-line. Over the years several theories have been developed for beam-like structures (see, for example, Love 1944, Antmann 1966, and Rubin 1997. This is because beam models have historically been used in many fields, from helicopter rotor blades in aerospace engineering to bridges components in civil engineering and surgical tools in medicine. Nevertheless, the interest in 30 advanced theories for complex beam-like structures has led to further researches also in recent years, due to the continuous need of ever more rigorous and application-oriented models. In this paper the attention is focused on the effects of important geometrical characteristics of those structures, such as the curvature of their center-line, as well as the twist and the taper of their cross-sections. After an introduction to modeling approaches for structures of this kind (section 2), a suitable model for the problem at hand is proposed (section 3). Finally, analytical results and numerical examples obtained by applying the 35 proposed modeling approach to reference beam-like structures are illustrated (sections 4 and 5).

Overview of modeling approaches
Aero-elastic modeling of modern blades can be addressed by means of different approaches (Wang 2016a). Those ones based on 3D FEM and beam models are two main choices for the structural part of this modeling. Although 3D FEM approaches can be very accurate and flexible, they can be computationally expensive for the analyses of complex systems, 40 especially if CFD analyses are executed in parallel. The overall computational cost can be reduced if faster aerodynamic models are used, such as the blade element momentum (BEM) model, but even this solution may not be efficient enough for aero-elastic analyses and multi-objective optimization tasks. The coupling of BEM and suitable beam models can be the best compromise between computational efficiency and accuracy. In this work we focus the attention on mathematical models to simulate the mechanical behaviour of complex beam-like structures (hereafter referred to as beam-like models, or BLM), 45 which can be curved, twisted and also tapered in their unstressed state, undergo large deflections, in-and out-of-plane crosssectional warping and small strain. Suitable models are needed to simulate the behaviour of those structures. In general, classical beam models (see, for example, Love 1944), which include extension, twist and bending, as well as the Reissner's formulation (1981), also accounting for transverse shear deformation, may not be sufficient. Geometrically exact models are a better choice, but a way to put them into a suitable form for engineering applications is usually needed (Antman 1966). In 50 general, suitable models should be both rigorous and application-oriented, two important requirements pursued over the years by many investigators (e.g. Giavotto 1983, Simo 1985, Ibrahimbegovic 1995, Ruta 2006, Pai 2011, and Yu 2012. For over a century researchers have sought to represent beam-like structures by means of 1D models. Several theories have been developed, from the elementary Euler-Bernoulli theory, to the classical theory which includes Saint-Venant torsion, up to more refined theories, such as the Timoshenko theory for transverse shear deformations, the Vlasov theory for torsional 55 warping restraint, and 3D beam theories which include 3D warping fields. Broadly speaking, beam theories can be grouped into engineering and mathematical theories. Several engineering theories are based on ad-hoc corrections to simpler theories (e.g. Rosen 1978), others are based on geometrically exact approaches (such as Wang 2016b and Hodges 2018). Among the mathematical theories, some approaches are based on the directed continuum (Rubin 1997), some others exploit asymptotic methods (Yu 2012). The reason for such a large amount of works on beam theory is that it has always found application in 60 many fields. By way of example, many approaches have been developed for helicopter rotor blades with an initial twist, but pre-twisted rods have always attracted the interest of many researchers in several fields (Rosen 1991). In the 1990's, Kunz (1994 provided an overview on modeling methods for rotating beams, illustrating how engineering theories for rotor blades evolved over the years. In those same years, Hodges (1990) reviewed the modeling approaches for composite rotor blades, discussing the importance of 3D warping and deformation coupling. More recently, Rafiee (2017) discussed vibrations 65 control issues in rotating beams, summarizing beam theories and complicating effects, such as non-uniform cross-sections, initial curvatures, twist and sweep. It seems that, unlike the case of the pre-twisted rods, published results for curved rotating beams with initial taper and sweep are quite scarce, although all these geometrical characteristics can play an important role.
Up to now much has been done to develop powerful beam theories. However, there is still a gap between existing theories and those that could be suitable for complex beam-like structures. In general, the geometry of the reference and current states 70 must be appropriately described. The curvature, twist and taper are important design features and should be explicitly included in the model. The analysis should not be restricted to small displacements. The model should provide the strain and stress fields in the three-dimensional beam-like structure, be rigorous and usable by engineers, and provide classical results when applied to prismatic isotropic homogeneous beams. Following these guidelines, a mathematical model to simulate the mechanical behavior of the considered beam-like structures is proposed hereafter. 75

Mechanical model for complex beam-like structures
Here we are concerned with developing a mathematical model to describe the mechanical behavior of beam-like structures which are curved, twisted and tapered in their reference state and undergo large displacements. One of the main issues with such a task is how to describe the motion of the structure (see, for example, Simo 1985, Ruta 2006, Pai 2014. The approach considered in this work is to imagine a beam-like structure as a collection of plane figures (i.e. the cross-sections) along a 80 regular and simple three-dimensional curve (i.e. the center-line). We assume that each point of each cross-section in the reference state moves to a position in the current state through a global rigid motion on which a local general motion is superimposed. In this manner, the cross-sectional deformation can be examined independently of the global motion of the center-line. So, it is possible to consider the global motion to be large, while the local motion and the strain may be small.

Kinematics and strain measures 85
We begin by introducing two local triads of orthogonal unit vectors. The first one is the local triad, b i , in the reference state, with b 1 aligned to the tangent vector of the reference center-line. This frame is a function of the reference arch-length s, that is b i =b i (s). The second local triad, a i , is a suitable image of the local triad b i in the current state. This frame is a function of the reference arch-length s and the time t, that is a i =a i (s,t). In general, a 1 is not required to be aligned to the tangent vector of the current center-line. See Figure 1. It shows a schematic representation of the reference (left) and current (right) states of a 90 beam-like structure. The generic cross-section in the reference state is contained in the plane of the vectors b 2 and b 3 . In the current state, if the cross-section remains plane (i.e. un-warped), it can belong to the plane of the vectors a 2 and a 3 . But the generic cross-section may not remain plane. So, we consider that its current (warped) state is reached by superimposing an additional motion to the positions of the points of the un-warped cross-section, as in Figure 1 (right).

Figure 1: Schematic of the reference and current states, center-lines, cross-sections and local frames
We continue by introducing the kinematical variables we use to describe the motion of the considered structure. To this aim, the orientation of the frames a i and b i relative to a fixed rectangular frame, c i , are defined as follows where A and B are two proper orthogonal tensor fields (i.e. their determinant is 1, see, for example, Gurtin 2003). We also 100 introduce a tensor field, T, which defines the relative orientation between the frames a i and b i as follows We define two skew tensor fields, K A and K B , and their axial vectors, k A and k B , which are related to the curvatures of the center-line of the structure in the current and reference states, as follows (see, for example, Simo 1985 andGurtin 2003) , , (3) 105 where the prime denotes derivative with respect to the arch-length, s, while the operator ∧ is the usual cross-product. In a similar manner, we introduce the skew tensor field Ω, and its axial vector field ω, related to the variation of the vectors a i over the time, t, as follows , where the operator ϕ[] provides the axial vector of the skew tensor between brackets.
The function R 0B , which maps the points of the center-line in the reference state, does not depend on time, while R 0A can change over the time t. Its variation is the time rate of change of the position of the points of the current center-line We are now in a position to introduce two important kinematic identities where γ and k are 00 It is worth nothing that γ and k vanish for rigid motions and are invariant under superposed rigid motion, hence, they are 120 well-defined measures of strain for beam-like structures (see, for example, Ruta 2006 andRubin 2000). Now, we start modeling the motion of the points of the cross-sections. In particular, we introduce two mapping functions, R A and R B , to identify the positions of the points of the 3D beam-like structure in its current and reference states. For what the reference state is concerned, we define the (reference) mapping function where R 0B is the position of the points of the reference center-line relative to the frame c i , b α are the vectors of the reference local frame in the plane of the reference cross-section, x α identify the position of the points in the reference cross-section relative to the reference center-line, and, finally, z i are independent mathematical variables which do not depend on time. In particular, z 1 is equal to the arch-length s, and z α belong to a bi-dimensional mathematical domain which is used to map the position of the points, x α , of the cross-sections. Throughout this paper, Greek indices assume values 2 and 3, Latin indices 130 assumes values 1, 2 and 3, and repeated indices are summed over their range.
It is worth noting that x k may or may not be equal to z k . The first choice leads to the most common modeling approaches available in the literature (see, for example, Simo 1985, Pai 2011, and Yu 2012. In this work we choose a set of relations between the position variables x k and the mathematical variables z k to provide a description of the shape of the considered beam-like structure, which can be curved, twisted and also tapered in its reference state. In particular, the span-wise variation 135 of the shape of the cross-sections is analytically modeled by means of a mapping of this kind where the coefficients Λ Bij are functions of z 1 . In the following we will consider the case of the curved and twisted beam-like structures with bi-tapered cross-sections, in which case the map (10) reduces to where the coefficients λ xα are functions of z 1 . It is worth noting that a suitable choice of those functions gives the possibility to reproduce interesting shapes. See, for example, Figure 2. It shows a 3D beam-like structure whose center-line is curved, while the cross-sections are twisted and tapered from the root to tip. The reference cross-sections in Figure 2 are ellipses with different sizes and orientations, but any other reference cross-section shape can be considered, such as the aerodynamic profiles which are commonly used for wind turbine blades, steam turbines blades, and helicopter rotor blades as well (see 145 also Griffith 2011, Bak 2013, Tanuma 2017, and Leishmain 2006 for examples of such profiles).

Figure 2: Example of curved, twisted, and tapered beam-like structure and local frame (left), taper and twist functions (right)
The position of the points in the current state are defined in a similar manner by means of the (current) mapping function where R 0A is a function mapping the position of the points of the center-line in the current state, while w k are the components of the 3D warping displacements in the local frame a k . The main formal difference between the reference and current maps is due to the warping, w, introduced to describe the geometry of the deformed state without a-priori approximation.
By using the maps (9) and (12), we can determine the 3D tensor, H, expressing the gradient of the current position, R A , with respect to the reference position, R B , as follows (see, for example, Rubin 2000) 155 When H is known, the 3D Green-Lagrange strain tensor, E, can be calculated (see, for example, Rubin 2000 andGurtin 2003). Hereafter we write the tensor E in a form based on simplifying assumptions applicable to the considered beam-like structure. In particular, we introduce the characteristic dimension of the cross-sections, herein denoted as h, the longitudinal dimension of the center-line, herein denoted as L, and we assume h to be much smaller than L. Moreover, we consider a thin 160 structure and assume the curvatures of its reference center-line are much smaller than 1/h (see also Rubin 2000). In addition, we assume the warping displacements, w k , are small. More precisely, by introducing a non-dimensional parameter ε much smaller than one, they are considered of the order of hε, while the order of magnitude of their derivative with respect to z 1 is εh/L. In general, all deformation measures, i.e. the 1D strain measures γ and k and the components of the 3D strain tensor, E, are assumed to be small. In particular, their order of magnitude is at most ε. For the considered structure, in the case of small 165 strains and small local rotations, we write the strain tensor, E, in the following form Let's now calculate the components of E by using (14) and neglecting terms smaller than ε. Algebraic manipulations only, too lengthy to be reported here in full, yields the following expressions for bi-tapered cross-sections In (15) Λ B22 and Λ B33 are the edge-wise and flap-wise taper coefficients (see, for example, Figure 2), while the components of the strain tensor, E, are taken with respect to the reference local frame, b i , i.e.
where • is the usual scalar (or dot) product and  is the tensor (or dyadic) product (see, for example, Rubin 2000).

Stress fields and constitutive models 175
Given the strain tensor, E, the stress fields in the structure can be calculated when a constitutive model is chosen. Limiting our attention to elastic bodies, in a pure mechanical theory, in the case of small strain, we use the following linear relation between the second Piola-Kirchhoff stress tensor, S, and the Green-Lagrange strain tensor (see, for example, Gurtin 2003) 2 where μ and λ are known material parameters related to the Young's modulus and Poisson's ratio.

8
In the case of small strains and small local rotations, we also write , T P TS C TST (18) where P is the first Piola-Kirchhoff stress tensor and C is the Cauchy stress tensor (Gurtin 2003). It is worth noting that in the considered case the tensor field T is sufficient to determine two of the above mentioned stress tensors (e.g. P and C) when the other one (e.g. S) is known. 185 We are now in the position to define the cross-sectional stress resultants, namely the force F and moment M. Using the first Piola-Kirchhoff stress tensor (Gurtin 2003), in the case of small warpings, small strains and small local rotations, we write 11 , where Σ is the domain corresponding to the cross-section on which the integration is performed and By combining Eqs. (15)-(19), the force and moment stress resultants can be related to the geometrical parameters of the structure and the 1D strain measures (8). However, such relations are actually known if we know the warping fields w k . An approach to obtain suitable warping fields is illustrated in section 3.4.

Expended power and balance equations
To complete the formulation, we conclude with considerations on the principle of expended power and the balance equations 195 for the considered structure. For hyper-elastic bodies (Gurtin 2003), we write the principle of expended power in the form where p are surface loads per unit reference surface (A), b are body loads per unit reference volume (V), Φ is the 3D energy density function of the body, and v is the time rate of change of the current position of its points, which is given by where the last term in (22) is the time rate of change of the warping displacement.
For small warpings, small strains, and small local rotations, if the power expended by surface and body loads on the warping velocities is neglected, the external power, Π e , reduces to the following form where the vector field v 0 is the time rate of change of the position of the points of the current center-line, the vector field ω is 205 the time rate of change of the orientation of the vectors a i , the terms F s and M s are suitable resultants of inertial actions and prescribed loads per unit length in the reference state, while the symbol Δ simply means that the function between brackets is evaluated at both the ends of the beam and the difference between those values is taken.
The cross-sectional warpings may be important in calculating the 3D energy function and cannot be neglected in the internal power, Π i . However, the internal power may be reduced to a useful form for beam-like structures by introducing a suitable 210 1D strain energy function, U. For example, if U can be expressed in terms of the strain measures, γ and k, we can write ( , , ) i ss where the vector fields f and m are defined in terms of the force and moment stress resultants, F and M, as follows ,

TT f T F m T M 
By using the principle of expended power, we also obtain balance equations for the vector fields F and M in the form 215 At this point, we have kinematic equations, (6)- (7), strain measures, (8) and (14), force and moment balance equations, (26), and the principle of expended power, Π e =Π i , in a suitable form for beam-like structures, (23)-(24). To complete the formulation of the model we need relations providing the 1D stress resultants in terms of the 1D strain measures. To this end, we need to know the warping fields. An approach to obtain suitable warping functions is discussed hereafter. 220

Warping displacements
In general, a 3D nonlinear elasticity problem can be formulated as a variational problem. In any case, if we try to solve the variational problem directly, the difficulties encountered in solving the elasticity problem remain. For beam-like structures whose transversal dimensions are much smaller than the longitudinal one, assumptions based on the shape of the structure and the smallness of the warping and strain fields can lead to useful simplifications. In particular, the resolution of the 3D 225 nonlinear elasticity problem can be reduced to the resolution of two main problems. See, for example, Berdichevsky (1981), who seems to be the first in the literature to plainly state that for elastic rods. One of those problems governs the local distortion of the cross-sections and is here referred to as the cross-sections problem. The other problem governs the global deformation of the center-line and is here referred to as the center-line problem. Hereafter, we consider the following variational statement to determine the warping fields which are responsible of the deformation of the cross-sections 230 In (27) the symbol δ stands for the variation operator and the density function Φ depends on the warping displacements. The warping fields satisfying (27) can be obtained by the corresponding Euler-Lagrange equations (see, for example, Courant 1953), by using numerical methods, in general, or analytical approaches providing closed-form expressions, in some particular cases. Once such a problem is solved, the components of the stress resultants (19) can be linearly related to the 235 components of the 1D strain measures, by using equations (14)-(19). Then, if it is preferred or deemed useful, those relations can also be arranged in a standard matrix form.
Note that to determine the current state of the structure we also need the displacements of its center-line points. They can be determined by solving the center-line problem, which is a non-linear problem governed by the set of kinematic, constitutive and balance equations introduced in section 3 (in particular, we are referring to the constitutive model in section 3.2, which 240 relates stress resultants and strain measures, and the balance equations for the stress resultants in section 3.3).
In the next sections we show some analytical solutions (section 4) and numerical results (section 5) that can be obtained by applying the proposed modeling approach to some reference beam-like structures.

First analytical results for bi-tapered cross-sections
In this section we consider the case of a beam-like structure with bi-tapered elliptical cross-sections. For this case we can 245 provide analytical solutions in terms of warping fields, while for generic shapes (e.g. the aerodynamic profiles used in wind turbine blades, steam turbines blades, and helicopter rotor blades as well) the problem corresponding to (27) can be solved by using numerical methods. However, this is not surprising, since even in the classical linear theory of prismatic beams analytical solutions are available for a limited number of cases only (see, for example, Love 1944).
As discussed in section 3, we are assuming the smallness of the warpings, strains and local rotations. Moreover, hereafter we 250 choose the current local frames to be tangent to the current center-line and we include possible shear deformations within the warping fields. In addition, with the aim of showing a first analytical solution for bi-tapered cross-sections, in this section we neglect the effects of the initial twist of the cross-sections. Then, we proceed to find a solution that satisfies (27).
Doing so, the Euler-Lagrange equations corresponding to (27), in the case we neglect the terms smaller than ε and maintain the terms related to extension, γ 1 , and change of curvatures, k i , are satisfied by the following warping fields 255 where d 2 and d 3 are the semi-major axes of a reference elliptical cross-section (e.g. the one at 18m from root section in Figure 2). Using this result, we can calculate the corresponding strain and stress fields, (14)-(18), stress resultants, (19), and strain energy function U. For example, if we consider a local frame in the reference cross-section with its origin at the crosssection's center of mass and its axes aligned with the cross-section's principal axes of inertia (as in Figure 2), we can write 260 the 1D strain energy function, U, in the form 2 2 2 4 2 3 4 2 4 2 1 1 1 2 2 3 3 In (29), E is the Young modulus, G is the shear modulus, while A, J 1 , J 2 and J 3 are the following geometrical parameters An interesting result is that the initial taper appears explicitly in all the previous relations (in terms of ρ and Λ). In its turn, 265 this allows an analytical evaluation of its effects on the 3D strain fields, which can be calculated by using (15) and (28), and which are required to determine the 3D stress fields (17).

Numerical simulations
In this section we provide the results of simulations conducted by using the modeling approach discussed in section 3, which we have implemented into a numerical code in MATLAB language. Those results are also compared with the results that can 270 be obtained by 3D-FEM simulations with the commercial software ANSYS.
In particular, we show a first set of test cases in which a beam-like structure with rectangular cross-sections undergoes large displacements, while it is fixed at one end and it is loaded at the other end by a force whose magnitude is progressively increased. In the second set of test cases we consider a more complex geometry, that is, a beam-like structure with elliptical cross-sections, which is curved, twisted and tapered in its reference configuration, while the loading condition is the same as 275 in the first set of test cases. Finally, in the third set of test cases, we consider (four) different beam-like structures under the same loading condition. In particular, we consider a first prismatic structure with elliptical cross-sections. The second structure is a modification of the first one, on which we maintain the same cross-section at 18m from root and we add the taper according to the taper coefficients of Figure 2. Starting from this latter, we consider a third structure which includes the twist of the cross-sections, assuming the twist law of Figure 2. The fourth one is a curved, twisted and tapered structure 280 obtained by the third one (tapered and twisted) by adding the center-line curvature. Then, we compare the results obtained by simulating the behavior of these four structures to shows the effects related to their geometrical differences.
In all the cases, the displacements of the points of the reference center-line are calculated by solving the center-line nonlinear problem through a numerical procedure we have implemented in MATLAB language, which is based on the kinematic, constitutive and balance equations introduced in section 3. In particular, we use the constitutive model introduced in section 285 3.2 to relate stress resultants and strain measures. We define the local frames orientation by using Euler angles and simulate orientation changes in terms of derivatives of those angles over the arch-length, s (see, for example, Pai 2003). We use the balance equations for the stress resultants introduced in section 3.3. Finally, we integrate (numerically) the resulting set of ordinary differential equations with respect to the arch-length, s. The results of this procedure are illustrated hereafter.

First set of test cases 290
In this set of test cases we consider a beam-like structure with rectangular cross-sections undergoing large displacements, while it is clamped at one end (i.e. the root) and it is loaded at the other end (i.e. the tip) by a force, F, whose magnitude is progressively increased (see Figure 3). The center-line length is d 1 =90m. The cross-section sizes are d 2 =8m (edge-wise) and d 3 =2m (flap-wise). The material properties are summarized by reference values of the Young's modulus, 70GPa, and 12 The simulations are run for different values of the tip force. The model we have implemented in MATLAB language to solve the non-linear problem provides results on the deformed configuration of the structure (e.g. Figure 3, left) within a few seconds. In all the cases, the simulation time is less than 2.4s. It is significantly less than that required by the corresponding non-linear 3D-FEM simulations carried out on the same computer, while the accuracy of the results is almost the same. A summary of the obtained results, in terms of global displacements and simulation times, is shown in Figures 3 and 4. 300 In particular, Figure 3 (left) shows the un-deformed shape (for F=0), as well as the deformed shapes for F equal to 10000kN, 25000kN and 50000kN. Figure 3 (right) shows the 3D-FEM deformed shape for F=25000kN (the coloured legend is for the flap-wise displacements). Then, Figure 4 (left) provides a comparison between the tip displacements obtained with the linear 3D-FEM, the nonlinear 3D-FEM and our model (therein referred to as 3D-BLM). It also shows the differences (between the non-linear 3D-FEM and the 3D-BLM) in terms of tip displacements and simulation times for the considered cases. 305

Second set of test cases
Let's now consider a more complex beam-like structure, more precisely, a 90m curved center-line with constant curvatures, which schematizes a pre-bent and swept beam whose tip is moved 4m edgewise and 3m flapwise, as in Figure 2. The local frames in the reference state are characterized by a pre-twist of 20deg/m. The reference cross-section at 18m from root is an ellipse whose semi-major axes are d 2 =2m (edge-wise) and d 3 =0.5m (flap-wise). The sizes of the other cross-sections change 315 according to the taper coefficients of Figure 2. For what the material properties are concerned, they are summarized by reference values of the Young's modulus, 70GPa, and Poisson's ratio, 0.25. Finally, the structure is clamped at its root and it is loaded by a flap-wise tip force, F, which varies from 100kN to 1000kN.
The simulations are run for different values of the tip force. The model we have implemented in MATLAB language to solve the non-linear problem provides the results about the deformed configurations (see, for example, Figure 5, left) within less 320 than 2.7 seconds. As in the first set of test cases, the simulation time is significantly less than that required by the nonlinear 3D-FEM simulations (they differ by at least one or two order of magnitude), while the accuracy of the results is again almost the same. Hereafter, we continue by showing some other information our model can provide. In particular, we can obtain the displacement fields of the points of the reference center-line ( Figure 6), as well as the change of curvatures of the beam-like structure (Figure 7, left) and the corresponding moment stress resultant (Figure 7, right). The moment components are in the 325 current local frame, a i , whose orientation has been determined as part of the solution of the nonlinear problem. For example, the orientation of the current local frame, a i , can be provided in terms of a set of Euler angles. See Figure 8. In this case we have considered the set of Euler angles corresponding to a first rotation, θ, about the initial z-axis, a second rotation, γ, about the intermediate y-axis, and a third rotation, ψ, about the final x-axis.

Third set of test cases
Here we consider different beam-like structures, starting with a prismatic elliptical one, including step by step the taper and 340 twist of the cross-sections and, finally, the curvature of the center-line, as discussed in the beginning of section 5. Note that the "curved-twisted-tapered" case considered here coincides with that discussed with more details in section 5.2 (see Figures   5-8, F=250kN). We proceed by simulating the bahavior of these four structures under a flap-wise tip force of 250kN. Then, we analyze the obtained results to show the effect of their geometric differences on their mechanical behaviour. A summary of the obtained results is hereafter. In particular, Figure 8 shows the reference and deformed states of the prismatic structure 345 (left) and the deformed states of the non-prismatic ones (right), while Figure 9 shows the displacements of their center-lines points. The main effect of the considered tip force is a displacement in the z-direction, in all the cases, with a displacement in the y-direction that we have only for the cases "tapered-twisted" and "tapered-twisted-curved", as it is expected. The obtained results have been compared with those of the 3D-FEM commercial software, for this third set of test cases too, confirming the computational efficiency and accuracy of the previous sets of test cases. A summary of those results is shown 355 in Figure 10, which provides a comparison in terms of tip displacements (components) for the four cases considered here.

360
As verified by many simulations and shown in the examples, the proposed approach performs well in terms of computational efficiency and accuracy. It can be used to study the mechanical behavior of beam-like structures, which are curved, twisted and tapered in their reference unstressed state and undergo large global displacements. It can provide information on the deformed configurations of those structures, such as their displacement fields, as well as the corresponding strain and stress measures. It is worth noting that it is suitable for beam-like structures with generic reference cross-sections shapes. However, 365 as already pointed out, for bi-tapered elliptical cross-sections we have analytical solutions in terms of warping fields, while for generic reference cross-sections shapes the problem (27) has to be solved by using numerical methods.

Conclusions
Wind turbine blades, as well as helicopter rotor blades, steam turbine blades and many other engineering structures, can be considered complex (non-prismatic) beam-like structures, with one dimension much larger than the other two and a shape 370 that is curved, twisted and also tapered already in the reference unstressed state. Their mechanical behavior can be simulated through suitable 3D beam models, which are computationally efficient, accurate and explicitly consider the main geometrical design features of those structures, the large deflections of their center-line, and the in-and out-of-plane warping of their cross-sections. In this work, curved, twisted and tapered beam-like structures have been modeled analytically. Their main geometrical characteristics have been explicitly included in the model. The warping displacement has been thought of as an 375 additional small motion superimposed to the global motion of the local frames. The resulting model is suitable to simulate large deflections of the center-line, large rotations of the local frames and small deformation of the cross-sections. The strain tensor has been calculated analytically in terms of the geometrical parameters of the considered structures, the 1D strain measures and the 3D warping fields. An approach based on an energy functional and a variational statement have been used to obtain suitable warping fields. The principle of expended power for curved, twisted and tapered beam-like structures has 380 been discussed, along with the balance equations for the force and moment stress resultants. Finally, analytical results and numerical examples, which include comparison with 3D FEM simulations, have been presented to show the effectiveness of the proposed modeling approach and the information it can provide.