Generalized Beam Theory for Thin-Walled Beams with Curvilinear Open Cross-Sections

The use of the Generalized Beam Theory (GBT) is extended to thin-walled beams with curvilinear cross-sections. After defining the kinematic features of the walls, where their curvature is consistently accounted for, the displacement of the points is assumed as linear combination of unknown amplitudes and pre-established trial functions. The latter, and specifically their in-plane components, are chosen as dynamic modes of a curved beam in the shape of the member cross-section. Moreover, the out-of-plane components come from the imposition of the Vlasov internal constraint of shear indeformable middle surface. For a case study of semi-annular cross-section, i.e., constant curvature, the modes are analytically evaluated and the procedure is implemented for two different load conditions. Outcomes are compared to those of a FEM model.


Introduction
Thin-walled beams are structural elements with three characteristic dimensions of different orders of magnitude: the thickness is small when compared to the dimensions of the cross-section, which in turn are small when compared to the beam length.
Starting from the second half of the 19th century, slender thin-walled members have found applications in civil and naval engineering as beams, columns, frame-works etc. Later they have been adopted to the needs of aeronautics and aerospace design, as well as, e.g., turbomachinery industry. They offer a number of advantages over the classical compact section beams. The most important ones are lightweight combined with strong rigidity and loads resistance, lower manufacturing costs due to reduced resource consumption and labor-saving design, lower transport and maintenance expenses etc.. Moreover, thanks to their specific layout they leave to the designer much more flexibility regarding the choice of the material and cross-section shape to meet any specific design requirements.
The slender thin-walled elements exhibit significantly different kinematic behavior when compared to solid section beams. In particular the profile warping occurs combined with elastic couplings involving bending and twisting; other local effects arise as well. Therefore, the recognized classical beam theories like Euler-Bernoulli and Timoshenko cannot be trustworthy to analyse thin-walled members.
The first successful attempt to model the kinematics of thin-walled beams was done by Vlasov [1]. He developed the Theory of Sectorial Area postulating the non-uniform torsion along the beam axis contributing to out-of-plane warping of the cross-section. This hypothesis was accompanied were used to describe the displacement field of the member cross-section. The admissible set of trial functions was determined requiring that the classic Vlasov's kinematic hypotheses of the linear theory (i.e., (a) transverse inextensibility and (b) unshearability) were satisfied also in the nonlinear sense. An illustrative example of C-lipped cross-section beam with uniform thickness subjected to a uniform vertical pressure load was presented. The accuracy of GBT analytical results was confirmed by the outcomes of FE simulations.
Studying the invoked above references and other GBT related literature one observes the vast majority of research deals with polygon-like cross sections while the cases of non-prismatic open or closed profiles are very scarce. This results from the fact that curved geometries make the problem to find the cross-section deformation modes quite complex. The exception might be tubular beam designs studied by Schardt in [14], who presented solutions to the problem of uniform bending of a closed cylinder specimen. The analysis was enhanced by discussion of solution accuracy when approximating the cylinder profile by a regular polygon shape.
In recent years, Silvestre [26] investigated the buckling behavior of circular hollow section (CHS) members. He postulated to extend the cross section analysis and include axisymmetric and torsional modes apart from the set of classical orthogonal shell-type deformations of the profile. The presented discussion concerned the influence of member length on the critical stress magnitude and corresponding buckling mode shape. Comparison of the analytical GBT results to FE simulations by shell elements revealed very good consistency with respect to both local and global structural behavior. Afterwards, the author extended his research to account for elliptical cylindrical shells and tubes under compression [27]. The problem of variable curvature along the cross-section mid-line was solved by introducing parametric functions of tangential angle and their Fourier series expansions to approximate the local arch radius. The proposed formulation proved to be very effective since only three deformation modes (one local-shell, one distortional and one global) were required to accurately evaluate the buckling behavior of elliptical members for a wide range of specimen slenderness. The structural behavior of circular cross-section members was studied also by Nedelcu [28]. The author considered conical shells under three different boundary conditions and subjected to buckling load. In the proposed approach the member deformed configuration was approximated by a combination of the predetermined shell-type deformation modes.
Thereafter, Luongo and Zulli [29] adopted a GBT framework to analyse the ovalization of a tubular cross-section beam when subjected to bending. In the kinematic analysis of the member cross section, distortional and bi-distortional strains were introduced in addition to the regular strain measures of rigid cross-section beams. Several cases of static loadings applied at the free end were investigated followed by large-amplitude free vibrations analysis. This research was continued to account for tubular beams with a soft core, possibly made of foam materials [30]. The findings showed the inclusion of a foam core could improve performance of the pipes due to the shifted structural instability limit load. In the subsequent paper [31], Zulli adopted this model to capture the double-layered pipe designs. New terms representing longitudinal slipping of the two layers were added to the kinematic description of the structure profile. Next, a homogenization procedure was applied to obtain the nonlinear, coupled, elastic response function of the beam-like member.
An efficient method to model beams with arbitrary but open cross section shapes within Generalized Beam Theory was presented by Gonçalves and Camotim [32]. The core of this approach was a modified two-stage cross sectional analysis procedure. Within the first step the curved mid-line of the cross-section was finely approximated by a series of straight segments constituting the polyline. The second stage consisted in selection just a few of natural nodes of this polyline to generate the independent DOFs for warping analysis while the intermediate nodes of the profile were treated like in classic procedure. By this approach it was possible to describe the cross-section geometry accurately, without generating an excessive number of deformation modes. The method was applied to several cross-sections with curved walls and it was shown it assured sufficiently accurate results even with a rough selection of cross-section DOFs. The method proved to be particularly efficient for polygonal sections with rounded corners.
The presented above discussion and review of thin-walled structures bibliography demonstrates that the problem of structural behavior of thin-walled beams with curvilinear open profiles has not been comprehensively studied within the GBT framework. The main purpose of present contribution is to fill this gap and to derive the analytical model of a generic thin-walled curvilinear cross-sections beam within the GBT modelling framework. Therefore the outline of the paper is as follows. In Section 2.1 the analysis of the profile kinematics is presented and expressions for strains in curvilinear reference system are given. The constitutive law is introduced in Section 2.2, and the procedure for setting trial in-plane and out-of-plane warping functions is explained in Section 2.3. Next, in Section 2.4, the governing equations are derived by means of the Hamilton's method. In Section 3, the numerical examples of semi-annular specimen in two loading conditions are presented. The analytical results are compared to the outcomes of the FE simulations. Finally, the paper closes in Section 4 with the concluding remarks.

The Beam Model
A thin-walled beam made of isotropic, homogeneous and linearly elastic material is considered. The specimen is initially straight and its cross-sections, which are constituted by curved webs and flanges, are assumed open, while the closed and multi-cell cases are left for future developments. Furthermore, it is assumed that the profile of the beam is constant spanwise with no taper nor pretwist.
First, linear strain-displacement relationship is written and then, following the hypothesis of plane stress in profile walls, the constitutive law for linear elastic material is used to calculate strain energy. Finally, Hamilton's principle is used and the equilibrium equations in terms of kinematic descriptors are obtained, imposing that the total potential energy attains a minimum in the class of compatible displacements. According to GBT framework, the displacements are considered as linear combinations of pre-established trial functions and kinematic descriptors which change along the axis line of the beam. The adopted cross-section modal shapes describe both rigid displacements, as in classical beam models, and changes in the shape of the profile characteristic for thin-walled members.

Geometry and Kinematics
With reference to Figure 1, the axis line of the beam is defined in correspondence of the cross-sections centroids. It is a straight line given as: where0 is the position of the origin, which is the centre of the left tip cross-section, and x is an abscissa which runs in the interval [0, l], where l is the length of the beam; the unit vectors (e x , e y , e z ) are orthogonal to each other and define the canonical (global) basis. The initial positionx of a generic point P located at the mid-surface of the wall is: is the vector which defines the mid-curveγ of the cross-section, having curvilinear coordinate s as parameter. Specifically about the mid-curveγ of the cross-section, the parameter s is chosen so that dy ds 2 + dz ds 2 = 1. Furthermore, it is convenient to set-up a local coordinate system along the curve and given by the tangent, normal and binormal unit vectors (e s , e n , e b ) where-with respect to global basis-the relations follow: and the curvature around e b is defined as: Through the beam deformation, the point P undergoes a displacement u, which can be expressed in terms of components on both the global and local bases, namely u x e x + u y e y + u z e z and ue s + ve n + we b . According to Equation (4) the two sets of displacements are related as: In order to take into account the thickness h of the walls, the webs and flanges are considered as singly curved thin shells, whose mid-surface coincides withγ ( Figure 2). The displacements of any generic point located off the mid-surface (n = 0) are obtained from the mid-surface counterpart ones using Kirchhoff's thin plate assumption.
respectively. It is worth noting how the second term in the right hand side of Equation (8) results from cross-section curvature about binormal axis [33]. As a consequence, any point out of the mid-surface performs displacement d which has the following components on (e s , e n , e b ): respectively. Substituting Equations (7) and (8) into the above definitions, the local basis displacement components are as follows: The infinitesimal strain components within the curved thin wall are defined as [33]: where ε x , ε s are longitudinal strain along x, s, respectively, and γ xs is the shear strain. After substituting Equation (10) in Equation (11), the strain components attain the following expressions: (12) where the pure membrane strain components, i.e., relevant to the mid-surface (n = 0) and indicated with superscript M, and the flexural strain components, i.e., proportional to n and indicated with superscript F, are: and Following the idea of the Kantorovitch semi-variational method the individual components of the mid-curveγ points displacement are expressed as a linear combination of K triplets of trial functions (U j (s), V j (s), W j (s)), j = 1, . . . , K, and unknown amplitude parameters a j (x): (15) where the prime indicates derivative with respect to x. The functions (U j (s), V j (s), W j (s)) represent the three components of the assumed j − th deformation field, defined on the section profile.
The use of Equation (15) in Equations (13) and (14) allows one to express the strain measures in terms of the kinematic descriptors a j (x), for j = 1, . . . , K: and

Constitutive Law
It is assumed that the beam is made of an isotropic, homogeneous and linearly elastic material. Moreover, the plane stress state hypothesis is adopted to a thin shell representing any individual beam wall. For thin-walled structures this assumption is realistic and provides reliable results. According to this formulation, the two-dimensional stresses are functions of solely two coordinates x and s, while all the transverse stresses are negligible: Thus, the corresponding 3D Hooke's law is simplified to: where E, G, ν are the material Young's and Kirchhoff's moduli, and Poisson ratio, respectively. Equation (19), when solved in terms of stresses, becomes: Substituting Equation (12) in Equation (20), the expressions of the individual stress components can be decomposed to the sum of membrane and flexural related terms, as well: Finally, the use of Equations (16) and (17) allows one to express the stress components in (21) in terms of the trial functions (U j (s), V j (s), W j (s)) and kinematic descriptors a j (x), j = 1, . . . , K. The final expressions are omitted for the sake of brevity.

Trial Functions and Vlasov's Constraints
The regular Vlasov's internal constraints specific for any open cross-section require that the profile mid-curveγ is inextensible and the middle surface of the thin-walled beam is shear indeformable. Thereby two relations follow: As a consequence, regarding Equation (16), every j-th set of assumed trial functions must satisfy the following conditions: In this study, the trial functions (U j (s), V j (s), W j (s)), j = 1, . . . , K are evaluated in two stages as proposed by Ranzi and Luongo in [21]. At first, the in-plane deformation field components (U j (s), V j (s)) are evaluated as the dynamical normal modes of the planar unconstrained inextensional (i.e., satisfying the first Vlasov's condition) frame having the shape similar to profile mid-curveγ. Thereby the considered frame is constituted by a monodimensional curved beam element or assembly of them in case of piece-wise regularity ofγ. Having found the in-plane deformation modes, at the next stage the third trial function W j (s) is determined. This represents the component relevant to the out-of-plane warping of the cross-section and it is consistently evaluated by integration of Equation (24). The emerging unknown integration constant is obtained by imposing the orthogonality condition of W j (s) to the uniform axial extension, i.e., γ W j (s)ds = 0. Therefore, the triple of thus deduced deformation fields (U j (s), V j (s), W j (s)) satisfy both Vlasov's conditions given by Equation (22).
In the above analysis a uniform mass per unit length is assigned to the referenced frame in order to evaluate the in-plane deformation components; the latter must also be consistent with Equation (23), i.e., the curved beam constitutingγ is axially indeformable. Furthermore, a suitable normalization condition on the normal modes is imposed.
It is worth noting that, due to the lack of external constraints, among the set of determined in-plane modes, three independent rigid motions of the cross-section in its plane are present as well. These are two translations in directions e y and e z , respectively, and rotation about e x . Moreover, consistently with the internal constraint in Equation (24), the out-of-plane components relevant to the two translations in directions e y and e z turn out to describe rotations of the cross-section about axes e z and e y , respectively, as in an Euler-Bernoulli beam bending. Thus, it should be pointed out that, with the present approach, it is precluded the possibility of introducing independent bending rotations as it would be done in Timoshenko beam models.
An additional trial function, which describes the uniform extension in direction e x is defined as well, and it is assigned to ordering j = 1. Therefore, the triplets (U j (s), V j (s), W j (s)) for j = 1, . . . , 4 describe rigid motions of the profile and are given as follows: j = 1: (0, 0, 1); j = 2: ( dy ds , − dz ds , −y); j = 3: ( dz ds , dy ds , −z); j = 4: (y dz ds − z dy ds , y dy ds + z dz ds , − (y dz ds − z dy ds )ds + c 4 ). The other triplets (j > 4) describe in-plane deformation of the cross-section and its corresponding kinematically consistent out-of-plane warping.
It is worth mentioning that the Vlasov's conditions (22) can be removed in the framework of GBT, as well: This is the case, for instance, of structures where shear-lag might occur, as shown in [22]; there, further trial functions are introduced, which account for both shear and extensional deformation.

Equilibrium Equations
Equilibrium equations in terms of generalized amplitudes a j (x) are obtained imposing the stationary condition to the sum of elastic and load potential energy (conservative external loads are assumed). In particular, the elastic potential energy has the following expression: where V is the total volume of the beam and the differential volume element being dV = dndsdx.
In the above the membrane components of strain in tangential direction s and shear mid-surface indeformability are omitted according to Vlasov's conditions Equation (22). After expressing the stress and strain components in Equation (25) in terms of the unknown amplitudes a j (x), by means of the relevant equations from Sections 2.1 and 2.2, the elastic potential energy becomes: where a(x) := (a 1 (x), . . . , a K (x)) T , the dot stands for scalar product, the prime indicates derivative with respect to x and the expressions of the elements of the individual matrices are: It is worth noting that A, B, C are symmetric matrices. The load potential energy is evaluated for a general case of a force per unit area acting in correspondence of the mid-surface. The generally oriented force is defined by the three components vector p(x, s) = p x e x + p y e y + p z e z that yields the energy: [p x u x + p y u y + p z u z ]dsdx (28) where the displacement components can be expressed in terms of a(x) as well, by means of Equations (6) and (15). The minimum of the total potential energy is attained imposing that δ(U e + U l ) = 0, ∀δa j (x), j = 1, . . . , K, where δ is the variation operator, which produces the following set of Euler-Lagrange equations: and boundary conditions in x = 0, l: where the elements of f and g are: Equations (29) and (30) represent a linear, ordinary differential boundary value problem with constant coefficients, governing the evolution of the functions a j (x) along the span of the beam. Even if it can be analytically solved using the well-known method of variation of constants, which provides the exact solution, numerical tools capable to solve BVP will be used to find the solution, for practical purposes. In spite of this, the outcomes will be referred to as analytical solution in the following part of the paper.
It is noteworthy that the case of branched cross-section can be addressed piecewise evaluating alongγ the integrals in Equations (27) and (31).

Definition of the Case-Studies
A thin-walled cantilever of length l = 1 m is considered for numerical analysis. The cross-section is semi-annular, with radius R = 0.1 m, uniform curvaturek = 1/R = 10 m −1 , and constant thickness of the wall h = 0.003 m (see Figure 3a). Defining the curvilinear abscissa s from the bottom end of the semi-annular, the curveγ turns out to be described by:

Trial Functions
With reference to Section 2.3, the in-plane trial functions (U j (s), V j (s)) are analytically evaluated as the in-plane normal modes of a semi-annular inextensible Euler-Bernoulli beam matching the shape of the cantilever cross-section (Figure 3a), without any external constraints. Differently, in case of cross-sections with general expression of the curvature as a function of s, numerical techniques might be called for. Exclusively for this step, the mass per unit length m and the bending stiffness EI are assumed uniform and, conventionally, their values taken as m = 1 kg/m and EI = 1 Nm 2 . The normal modes of the arch beam turn out to be the solutions of the following eigenvalue problem (see [34]): with natural boundary conditions in s = 0, πR: where the above represent the normal and shear forces and bending moment, respectively. In particular, Equation (33) (36) is the classical constitutive law for bending moment M and curvature k b , and Equations (37)-(39) are the balance equations for normal N and shear T forces, and bending moment M, respectively. The eigenvalue λ is the square of the circular frequency ω. It is worth noting that no constitutive law is used for N and T, they being the internal reactive forces to the constraints (33) and (34).
Assuming an exponential solution of type v(s) =ve Ωs where v(s) = (U, V, Θ b , k b , N, T, M) T , after usual steps in modal analysis, the following bi-cubic dispersion equation is obtained: Its roots Ω i (λ), i = 1, . . . , 6, after substitution in Equations (33)-(39), allow one to evaluate the relevant eigenvectorsv i and build the solution v(s) = ∑ i c ivi e Ω i s , where c i are six undetermined constants. The latter expression is substituted in Equation (40) to get the first K natural frequencies (or λ) and, for each of them, five of the six coefficients c i , while the sixth is given by the normalization condition.
As expected, λ = 0 is a solution of multiplicity 3, giving rise to three in-plane independent rigid motions. An additional rigid motion, required to describe the uniform displacement in direction e x , is prepended to them as reported at the end of Section 2.3.
Once the modal in-plane components U j (s), V j (s), j = 1, . . . , K are obtained, the out-of-plane ones W j (s) are evaluated as well, as described in Section 2.3.
The first four modes, relevant to rigid motions, are shown in Figure 4 where, in the first and second row, the in-plane and out-of-plane components are drawn, respectively. In particular, in the first line, the in-plane components (U j (s), V j (s)) are combined to show the actual displacement of the cross-section (in thick solid line), whereas, in the second row, the corresponding normal direction component W j (s) is plotted (as polar plot, still in thick solid line). Following the same rule, the first four in-plane deformation modes (i.e., which correspond to λ = 0) are shown in the first row of Figure 5 and, in the second row, the corresponding out-of-plane components.

Equilibrium Configuration
Using the K = 8 trial functions discussed in Section 3.2, the equilibrium Equation (29) can be evaluated through Equation (27), and combined to the boundary conditions (30) relevant to the cantilever, i.e., geometric at x = 0 and natural at x = l. The obtained equations are then solved using a finite difference code based on the three-stage Lobatto IIIa formula, as implemented in the bvp4c library command in MATLAB [35], and the solution in terms of a j (x), j = 1, . . . K, is used to calculate displacement and stress components in specific points of the thin-walled beam.
In particular, for the two load conditions, plots are shown for: kinematic descriptors a j (x) as functions of x ∈ [0, l]; displacement components u x , u y , u z as functions of x for three points of the mid-curveγ, i.e., at positions where s = 0 (referred to as position S, i.e., the most southerly point), s = πR 2 (referred to as position W, i.e., the most westerly point), s = πR (referred to as position N, i.e., at the most northern point)-see also Figure A1  For the load condition 1 (load in e z direction), the solution is shown in Figure 6. From Figure 6a,b it can be observed that the component a 1 , corresponding to uniform extension, is always zero, as expected in a linear problem with structure subjected to transverse loads only. Moreover, amplitudes a j , j = 3, . . . , 6 have at least one order of magnitude larger than the remaining a 2 , a 7 , a 8 , revealing that global displacement along e y is almost zero and higher modes are actually not required for this case (it is verified, even if not reported here for the sake of brevity, that higher modes than K = 8 give vanishing contributions for the two analyzed cases). In particular, the components a 3 , a 4 , a 5 are those which give larger contributions, indicating that besides a rigid motion of the cross-section, a not-negligible deformation component related to a 5 (and in minor intensity a 6 ) is present.
In Figure 6c-e, the individual displacement components at the three points of the cross-section (S,W,N) are shown: the component u x comes only from warping of the cross-section and is much smaller than u z and u y . Significant in-plane rotation (about e x ) and change in shape of the cross-section are observed, as related to quite large a 4 and a 5 , respectively, as well as the expected slightly larger maximum displacement in direction e z (direction of the load). These outcomes reveal very good agreement with the FEM results. Consistently with Vlasov's theory of non-uniform torsion [1], Figure 6f   x (m) (e) For the load condition 2 (load in e y direction), due to the symmetry of the problem about axis e y , only modes respecting the same symmetry condition are used in the procedure, namely: II, V, VII in Figures 4 and 5. The relevant amplitude components a 2 , a 5 , a 7 are shown in Figure 7a with clearly dominating role of the component a 2 corresponding to rigid motion of the cross-section. Furthermore, among them, a 7 is almost negligible as well. In Figure 7b,c, the symmetry of the problem causes the blue and green lines to be superimposed to each other, and satisfying agreement with FEM outcomes is noted. Distortion of the cross-section is evident, as confirmed by the displacement components u z in Figure 7d, which are opposite in sign for points S and N, respectively, as given by a 5 multiplying the V-ip trial function in Figure 5. Specifically for this displacement component, a non vanishing difference of GBT analytical and FEM numerical results is detected at the free tip relevant to positions N and S. The outcomes are anyway in very good agreement in most of the domain of the beam, i.e., for x ∈ [0, 0.8]. Figure 7e shows the axial stress component σ M x . The distribution is fully consistent with the contribution expected in simple bending test, with stress proportional to a 2 and multiplying the component II-op in Figure 4. Finally, a vanishing change in the value of resultant σ x along the thickness of the cross-section is shown at positions S and N in Figure 7f. This is due to the fact that the thickness of the profile wall in these points is perpendicular to the load direction thus the flexural component σ F x is immaterial.

Conclusions
A thin-walled beam with open and curved cross-section is considered in the paper, with the aim of formulating an analytical Generalized Beam Theory framework model for this specific case. Consistently, the displacements of the points of the walls are expressed as linear combinations of trial functions and unknown amplitudes. Then, by adopting the Hamilton's principle, the equilibrium equations under external loads are evaluated taking into account the initial curvature of the members. As proposed in the literature [21], two components of trial functions are chosen as dynamical in-plane modes of the cross-section considered as a free beam. Consequently, the third out-of-plane component is evaluated imposing the validity of Vlasov's internal constraint. Being the member cross-section in the shape of an arch, the trial functions are analytically evaluated solving the relevant eigenvalue problem, which provides both rigid (frequency zero) and deformation (frequency different than zero) motions. In case of elastic, isotropic and homogeneous material a case-study subjected to two load conditions is considered. The outcomes in terms of displacements are compared to those of a companion FEM model, showing generally good agreement. The provided results highlight the crucial contribution of the cross-section change in shape at equilibrium.

Funding:
The research was financed in the framework of the project Regional Excellence Initiative funded by the Polish Ministry of Science and Higher Education (contract no. 030/RID/2018/19) and carried out at Lublin University of Technology.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. FE Model
To verify the analytical results, the ABAQUS system finite element representation of the structure was developed. The beam under consideration was modelled by conventional linear quadrilateral shape shell elements with reduced integration (S4R). They are general-purpose elements capable to account for shear effects and provide correct results of displacements, strains, and stresses [36]. This type of elements is recommended for both thick and thin wall specimens approximating the Kirchhoff model as the thickness decreases. Moreover, this element accounts for finite membrane strains and arbitrary large rotations; therefore, it is suitable also for geometrically nonlinear analysis.
Alternatively, one can use continuum shell elements but they require the geometry to model explicitly the thickness of the shell, as it would be done with classical 3D solid elements. However, in contrast to the classical shell, continuum shell elements enforce the first order shear deformation constrains through particular element interpolation functions.
The beam structure was meshed into 100 elements along the span and 30 ones along the profile perimeter. Moreover, the nodal sets required to record the spanwise deflections at most southerly (S), most westerly (W) and most northern (N) positions of the profile were defined as shown in Figure A1.