Isogeometric analysis of free-form Timoshenko curved beams including the nonlinear effects of large deformations

Seyed Farhad Hosseini·Ali Hashemian·Behnam Moetakef-Imani·Saied Hadidimoud

Abstract In the present paper,the isogeometric analysis(IGA)of free-form planar curved beams is formulated based on the nonlinear Timoshenko beam theory to investigate the large deformation of beams with variable curvature.Based on the isoparametric concept,the shape functions of the field variables(displacement and rotation)in a finite element analysis are considered to be the same as the non-uniform rational basis spline(NURBS)basis functions defining the geometry.The validity of the presented formulation is tested in five case studies covering a wide range of engineering curved structures including from straight and constant curvature to variable curvature beams.The nonlinear deformation results obtained by the presented method are compared to well-established benchmark examples and also compared to the results of linear and nonlinear finite element analyses.As the nonlinear load-deflection behavior of Timoshenko beams is the main topic of this article,the results strongly show the applicability of the IGA method to the large deformation analysis of free-form curved beams.Finally,it is interesting to notice that,until very recently,the large deformations analysis of free-form Timoshenko curved beams has not been considered in IGA by researchers.

Keywords Curved beams·Nonlinear Timoshenko beam theory ·Large deformation ·Isogeometric analysis·NURBS curves

1 Introduction

A beam is typically defined as a structural member with a large ratio of length to cross-sectional dimensions,which undergoes stretching along its length and bending about an axis transverse to the length[1].With a wide variety of applications in aerospace,civil and mechanical structures,beams are one of the well-known elements in most engineering analyses.Moreover,based on the geometrical aspects of modern design,curved beams with variable curvature are increasingly utilized by engineers and architectures in recent decades.

From the technical point of view,when the applied transverse load on the beam(straight or curved)becomes large,the linear load-deflection relationship is no longer valid,since the beam develops internal reaction forces that resist deformation.The magnitude of these internal forces increases with the loading and deformation[1,2].This nonlinear load-deflection behavior of beams is the main topic of this article.The nonlinearity in the formulation comes from the inclusion of second order terms in the strain-displacement relations.Asa result,the stiffness matrix in the finite element analysis(FEA)will itself be a function of the displacement field,in other words,[K(Δ)]{Δ}={F}where K, Δ,and F are the stiffness matrix,displacement field,and applied force,respectively.

Two different theories for modeling the kinematic behavior of beams are generally classified in the literature as the Euler–Bernoulli beam theory(EBT)and the Timoshenko beam theory(TBT).The main difference between these theories is in accounting for transverse shear strain,which is neglected in the Euler–Bernoulli theory.The classical beam theory is based on the Euler–Bernoulli hypothesis that plane sections perpendicular to the midline of a beam before deformation remain plane,rigid(do not deform),and rotate in such a way that they remain perpendicular to the midline after deformation.This assumption neglects the Poisson effect and transverse strains[1].Hence,in order for a beam to be analyzed more accurately,the Timoshenko beam theory should be applied.

In the structural analysis of curved beams,due to their complex geometrical forms,the deformation of the beam depends on the coupled equations between the tangential and transverse displacements,as well as rotation by curvature effects[3,4].As a result,analytical approaches can rarely find a closed-form solution for the displacement fields of curved beams,and a finite element analysis should be performed.However,this requires converting the geometries created by computer-aided design(CAD)into FEA models,which usually reveals certain shortcomings.This is because conventional FEA approximates the modelby linear or quadratic Lagrangian shape functions[5]while free-form curved beams are mainly represented by parametric nonuniform rational basic splines(NURBS)in CAD software[6].To address this issue,structural analyses are needed to analyze accurately these structures based on the NURBS formulation.One of very first researches in the literature concerning the incorporation of NURBS into FEA is conducted by Ganapathi et al.[7],who proposed a curved cubic B-spline beam element for static analysis.Among the other research,Hashemian and Imani[8]recently employed the B-spline representation to assess the automotive bodies whose deformations are computed by nonlinear FEA.However,the most comprehensive and highly-cited approach is isogeometric analysis(IGA),introduced by Hughes et al.[9],that has been developed to formulate finite elements for arbitrary free-form geometries.Isogeometric analysis is based on the concept of isoparametric formulation in FEA where the parametric representation of both the geometry and the displacement fields is of the same order(in contrast to super and sub-parametric formulations where the orders of shape functions of geometry and displacement fields are different).The IGA approach employs the same NURBS shape functions to represent both the geometry and unknown displacement fields in the FEA.Therefore,a seam less transition will be made between the CAD models and the structural analysis,while also eliminating the need for generating new meshes in the FEA[10–12].It is also interesting to note that the isogeometric approach has been applied to a wide range of engineering problems such as solid mechanics[13–18],composites and functionally graded materials[19,20], fluid–structure interaction[21],heat transfer[22],and eigenvalue,boundary value and initial value problems[23].

Recently,isogeometric analysis has been widely used in the structural analysis of curved beams.Bouclier et al.[24]investigated the use of parametric NURBS curves with higher degrees in order for the Timoshenko beams to be insensitive to shear and membrane locking issues.Locking phenomena commonly take place in FEA with shape functions of lower degrees,where the element cannot satisfy the strain constraints[1].The work of Cazzani et al.[25]also stated that membrane and shear locking phenomena can be controlled by either properly choosing the number of elements or the degree of the NURBS functions.Furthermore,Hosseini et al.[26]studied the effect of the parameterization of the NURBS shape functions on the isogeometric analysis of free-form curved beams.In their study,Luu,et al.[27]implemented the IGA in a vibration analysis of Timoshenko curved beams with variable curvature.However,all the aforementioned research is based on the assumption of a linear force-deflection relationship.

With special attention to large deformations of curved beams,Bauer et al.[28]derived the element formulation by means of nonlinear kinematics and utilized IGA to perform a finite element analysis.But their research was limited to the Euler–Bernoulli beam theory.It is interesting to notice that,until very recently,the isogeometric analysis of free form curved beams based on Timoshenko beam theory and including nonlinear effects of large deformations has not been considered by researchers.Hence,drawing attention to the central importance of the mentioned topic in practical engineering problems,the objective of this article is to improve the current IGA approach to take into account the nonlinear force-displacement relationship of free-form curved beams based on the Timoshenko beam theory.For this purpose,in Sect.2,the nonlinear Timoshenko curved beam theory is developed and compared to the straight beam theory.Then,starting with a brief introduction to NURBS curves in Sect.3,the isogeometric analysis is explained and the finite element formulation is derived.Section4 illustrates the applicability of the proposed approach by means of some case studies and gives a complete discussion of the results,followed by Sect.5,which draws some conclusions.

2 Non linear Timoshenko curved beam theory

The Timoshenko beam theory(TBT)is based on the fact that a straight line transverse to the axis of the beam would not remain normal to the midline after deformation since the rotation due to flexure is independent of the slope of the beam[1].In this section,the nonlinear Timoshenko beam theory isfirst described for straight beams and then extended to curved to curved beams with variable curvature.For both straight and curved beams,displacement fields are assumed for the beam and strains that are consistent with the kinematic assumptions of the Timoshenko theory are computed.Then,a weak formula-tion is developed using the principle of virtual displacements.As stated earlier, including large deformations in the analysis makes the strain-displacement relations nonlinear.

Fig.1 Kinematics of a straight beam based on Timoshenko beam theory

2.1 Straight beam theory

The displacement fields of a straight Timoshenko beam,as depicted in Fig.1,can be given as follows[1]:

The nonlinear strain-displacement relations including the second order terms are,therefore,defined as follows[1]where γ is the transverse shear strain neglected in the Euler–Bernoulli beam theory.The first three terms in the equation ofε express the membrane strain while the last term accounts for the bending strain of the beam.

Fig.2 Kinematics of a curved beam based on Timoshenko beam theory

In the above equations,the termarising from the longitudinal stretching of the beam is often neglected for straight beams[1].However,negligence of this term may cause simulation error,especially in the curved beams when the quadratic term of the gradient of the longitudinal deformation approaches the quadratic term of the gradient of the transverse deformation[4,29].

2.2 Curved beam theory

Generalizing the kinematics of straight beams,the displacement field of a curved beam in the curvilinear coordinate system(s-z)is demonstrated in Fig.2 where ρ is the radius of curvature,u and ware tangential and normal displacements of the midline,respectively,and w′=d w/d s.

It is suggested by Babu and Prathap[30]to use the rotation due to flexure(ψ)as a nodal degree of freedom(DOF)for curved beam elements(similar to straight beams).In this case,the curvature and transverse shear strains are defined as,respectively.However,it is explained by Day and Potts[31]that ψ is not equal to the physical rotation of the section and,hence,is not suitable to be used as a nodal DOF for curved beams.The appropriate choice of nodal DOF,instead of ψ,would be ? = ψ +u/ρ[32].Consequently,the displacement fields of an arbitrary point on a free-form curved beam in the curvilinear coordinate system(Fig.3)are defined as

Fig.3 DOFs of a point on the midline of a free-form curved beam

By separating the membrane and bending(or curvature)strains and defining the curvature of the beam as κ=1/ρ,one can obtain the nonlinear strain-displacement relations as follows.Note that tow rite the equations more compactly,the superscript′refers to differentiation with respect to s.

2.3 Weak formulation

The weak formulation of structural problems can be derived using the principle of virtual displacements.This principle states that for a body in equilibrium,the virtual work done by external forces should be the same as the change in the elastic strain energy.In other words,

where E and V,defined by Eqs.(6)and(7),are the elastic strain energy and the work of external forces,respectively.

In the above equations,fs,fz,and m are the distributed longitudinal,transverse,and moment loads,respectively(Fig.4). Also,E is the Young modulus,G is the shear modulus,A is the cross section,I is the second moment of inertia,and l is the arc length of the beam.Finally,Ksis the shear correction coefficient.

Fig.4 External loads on a curved beam

Substituting Eqs.(4)into(6),the variational form of Eqs.(6)and(7)would be:

Fig.5 A cubic B-spline curve with seven control points

The weak formulation of the displacement fields of the curved beam can now be stated as follows,where Qs,Qz,and Qmare boundary loads.

3 Isogeometric analysis

3.1 NURBS curves:definition

NURBS curves play a central role in the isogeometric analysis of curved beams since they serve both for representing the geometry and for expressing the displacement fields.Starting with the definition of the curves,an NURBS curve of degree p,is a piecewise continuous function with n+1control points Pi=[Xi,Yi]and corresponding weights wi,expressed as follows[6]:

A simple form of an NURBS curve is a B-spline curve in which all weights are equal to 1. In this case, the summation of the basis functions in the denominator of Eq.(13)is also equal to 1 and,therefore,Ri,p(ξ)=Ni,p(ξ)or in other words,

Fig.6 Cubic basis functions of B-spline curve of Fig.5

For example,Fig.5 depicts a cubic B-spline curve with 7 control points and 4 non-zero knot spans,i.e.Ξ ={0,0,0,0,0.25,0.5,0.75,1,1,1,1}.This curve can be thought of as the geometry of a planar curved beam.As indicated in the figure,the first and last control points are coincident with the starting and ending points on the curve.The cubic basis functions,which are of order 2 continuous(C2)at knots,are demonstrated in Fig.6.

3.2 NURBS curves:geometry construction

It is interesting to note that if the data points representing the geometry of a curved beam are available,the NURBS expression of the beam can be found by curve fitting technique[6,33].Generally,given a set of r+1 data points Dk=?Dkx,Dky?where k=0,1,...,r,the geometry is to be constructed in such a way that the control points are the output of a global curve fitting problem.For this purpose,the parametercorresponds to k-th data point and should be determined by a parameter selection method uniformly spaced,chord-length or centripetal as represented by Eqs.(16)to(18),respectively:

where

and

In this case,the number of control points is to be selected in the way that the desirable fitting error will be achieved.It should be noted that NURBS curves approximation with unknown weights and control points leads to a nonlinear objective function[33,34]However,by selecting a suitable number of control points, it is possible to achieve good results in approximating by B-splines.Afterwards,the weights can be modified in order to iteratively change the shape of the curve if needed.It is worth mentioning that the internal knots in Eq.(20),to guarantee that every knot span contains at least oneshould be defined as follows wherei=int(jd)and α =jd ? 1[6]:

Figure 7 demonstrates some data points that are approximated by a cubic B-spline curve with 10 control points based on different parameter selection schemes.The figure shows that all curves obtained by these three parameterization strategies are almost geometrically the same.However,since they are parametrically different,they may amount to different convergence rates in the IGA solution(see Sect.4.6).

3.3 Isogeometric model of curved beams

In isogeometric analysis,both the geometry and displacement fields in FEA are represented by means of B-spline or NURBS basis functions.This means that the shape functions of the unknown field variables of interest(u,w,and ?)are considered to be the same as the basis functions defining the geometry.Table 1 presents the geometry and displacement fields of curved beams in IGA where x and y state the position of an arbitrary point on the midline of the beam,[Xi,Yi]is the position of the i-th control point in the x-y plane,Ui,Wi,andΦiare the control points of the field variables,and Ri(ξ)is the i-th bas is function of degree p as introduced in Sect.3.1.

Fig.7 Geometry construction from data points by a cubic B-spline curve with 10 control points based on different parameter selection schemes

Table 1 NURBS representation of the geometry and displacement fields of curved beams

It should be pointed out that the knots on the NURBS curve(Fig.5)are equivalent to the positions of finite element nodes in IGA.Hence,an element in the FEA corresponds to a non-zero knot span in the parameter space of the NURBS geometry.This means that by increasing the number of control points,or the degree of the curve,or both,it is possible to refine the mesh size on the beam.These refinements are referred to as h-,p-,and k-refinements,respectively[27].

In isogeometric analysis,the transformation from the curvilinear coordinate s on the beam geometry to the curve parameter ξ,as schematically shown in Fig.8,is performed by the Jacobian,defined as

The curvature of the curved geometry can,therefore,be given as:

noting that by rewriting Eq.(13)as C(ξ)=A(ξ)/W(ξ),the k-th derivative of x and y with respect to ξ in the NURBS notation can be extracted from the k-th derivative of C(ξ)as follows[6]:

Substituting the displacement fields of Table 1 in the weak forms of Eqs.(10)–(12)and integrating over the element domain,the force–displacement relation[K]{Δ}={F}w ill be obtained:the element level are as follows,whereand leis the length of an element:

The explicit forms of the stiffness matrix and load vector at

In the above equations, the integration variable s can be converted to the parameter ξ using Eq.(23)as d s=J dξ.Hence,each term of the stiffness matrix and load vector can be computed using the standard Gauss–Legendre quadrature as a numerical integration technique.For this purpose,another transformation from the parameter space to the parent element is needed as expressed by Eq.(30)and depicted in Fig.8.In this paper, the selective and reduced numerical integration technique[35]is utilized to calculate the integrals.It is important to mention that quadratic B-spline curved beam elements may represent shear or membrane locking issues.However,by employing higher order basis functions for geometry and unknown field variables,the locking problems are generally alleviated[27,36].

Fig.8 Transformation from the geometric domain of the beam to parametric domain and parent element

Having included large deformations in the analysis of Timoshenko curved beams,the components of the stiffness matrix of Eq.(28)are functions of the unknown displacement field variables u and w.In other words,the force–displacement relation of Eq.(27)becomes[K(Δ)]{Δ}={F}.Therefore,an iterative method for solving nonlinear system of equations should be employed to find the unknown displacements.

3.4 Iterative solution of nonlinear equations

Generally,the system of nonlinear equations can be linearized using either the direct iteration or the New ton–Raphson iterative method.In this paper,the New ton–Raphson method is preferable since it converges faster.In this method,the solution at the r-th iteration is determined from the assembled set of linearized equations as follows where{R}is the residual vector and[T]is the tangent stiffness matrix[1]:

In this equation,the stiffness matrixis evaluated at the element level using the known solutionat the(r?1)-th iteration.The components of the tangent stiffness matrix?Te?at the element level are calculated using the following definition[1].Note that referring to Eq.(27),α,β =1,2,3,Δ1=u,Δ2=w,Δ3= ?,and the range of k is dictated by the size of the matrix?Kαβ?.

Finally,values of u,w,and ? are determined after the residual reaches the predefined value of error.

4 Case studies

In this section, five examples of various geometries and loading scenarios are presented.In the first example,a straight beam under uniformly distributed loading is presented.A convergence study is also performed for this example and the difference between the Timoshenko and Euler–Bernoulli beam theories are highlighted in results.The shallow-arch beam under amiddle point load is the second example, which represents a curved beam with constant curvature.A Tschirnhausen beam under constant distributed load that represents a free-form curved beam with variable curvature is investigated as the third case study.Finally,the cardioid beam under a point load and a logarithmic spiral beam under a bending moment,as two more free-form structures,are investigated.It should be pointed out that in the free-form examples,the NURBS geometries for the isogeometric analysis are constructed using the curve fitting technique described in Sect.3.2.Additionally,the effect of parameterization strategies in creating suitable geometries for IGA is investigated for these examples.For the first two case studies,which are directly modeled by the NURBS definition,the degree elevation and knot insertion algorithms are employed to create the required geometries as the mentioned algorithms do not change the shape of the curves[6].

4.1 Straight beam under uniform distributed loading

A clamped–clamped straight beam under uniform distributed loading is considered as the first case study of this article(Fig.9).It is important to note that the shear deformation effect in thick beams is the main reason for the difference between Timoshenko and Euler–Bernoulli theories.However,when the beam is thin,the Euler–Bernoulli approximation is sufficient and both theories give almost the same results.It will be discussed in detail at the end of this subsection.

Fig.9 A clamped–clamped straight beam under uniform distributed loading

In order to assess the validity of the presented IGA formulation,we start with abeam of length L=100 in(1 in=2.54 cm)and 1 in×1 in cross section.The beam is made of steel with elastic modulus and Poisson’s ratio of E=30M psi(1 M psi=6.895 ×103MPa)and ν=0.25,respectively,and subjected to a uniformly distributed load of magnitude q0.The maximum deflection at the center of the beam for different values q0is reported in Table 2.A cubic NURBS geometry with different numbers of control points is used for isogeometric modeling of this example.The convergence is reached after a maximum of five iterations in the most critical case(maximum nonlinearity).The results of this classical example are also compared to those presented by Reddy[1].

The table shows that when the number of control points is not large enough,the accuracy of the IGA solution decreases as the magnitude of the load is increased.However,for larger number of control points(i.e.,for more number of elements in IGA),the error is significantly decreased that shows the accuracy of the IGA approach.This is because the nonlinearity increases with the applying force and a model with the lower number of control points(i.e.,the lower number of elements)cannot express such nonlinearity very well.Figure 10 also draws the same conclusion for different numbers of control points.The comparison of linear and nonlinear load-displacement results are also shown in the figure.The results of IGA closely follow the nonlinear FEA results of Reddy[1],which proves the validity of the presented method of this article.

The difference between Timoshenko and Euler–Bernoulli beams can be highlighted by considering the effect of the length to thickness ratio(L/H)on the non-dimensional deflection.The non-dimensional deflection for the unit width of the beam is defined as follows[1]:

where wmaxis the maximum deflection at the center of the beam.Figure 11 illustrates the load-displacement response predicted by IGA for TBT in the current study and compares the results to those of EBT.In this figure,the shear deformation effect in the Timoshenko theory is the main reason for the difference in the results when the beam is thick(L/H=10).However,when the beam is thin(L/H=100)the Euler–Bernoulli approximation is sufficient and both theories give almost the same results.

4.2 Large deflection analysis of shallow-arch beam

The large deflection of a clamped circular arch with a concentrated load at the middle is analyzed as the second example.All dimensional and material data are provided in Fig.12.For this example,one half of the arch is modeled due to symmetry.

Table2 Maximum transverse beam deflection(in inches)for first case study:load values versus number of control points

Fig.10 Comparison of linear and nonlinear load-displacement results for different numbers of control points in IGA

Fig.11 Effect of length to thickness ratio(L/H):Timoshenko versus Euler–Bernoulli beam theories

Fig.12 A clamped–clamped shallow-arch beam under concentrated load

Fig.13 Comparison of large deflection analyses of shallow-arch beam with present formulation of this article

The vertical displacement at the apex(point of load)is considered as the output of the analysis.Since a symmetrical modeling is adopted throughout this example,the apex is coincident with the ending point of the NURBS curve representing the geometry of one-half of the beam.It was stated in Sect.3.1 that the ending point of the curve is equal to the last control point.Therefore,the point load can be imposed exactly on the last control point as a boundary load.

The results of this well-established benchmark example are superimposed onto a diagram created by Bathe and Bolourchi[37]and demonstrated in Fig.13.Four isogeometric analyses are performed with different degrees and numbers of control points and compared to the results of Mallett and Berke[38]as a reliable solution to this example.The results of Dupuis et al.[39],who analyzed the same arch using curved beam elements using their“Lagrangian”and “updated”formulations,are also illustrated in the figure.This comparison proves the validity of the developed model in nonlinear analysis of curved beams with constant curvature.It should be noted that the applied load could be increased up to around 35 pounds after which the buckling phenomenon would take place[40].However,our analysis is performed for nonlinearities in the pre-buckling region and the main purpose is to compare the results with abovementioned research works.

Fig.14 Clamped–clamped Tschirnhausen beam under uniform distributed loading

Fig.15 Convergence rates of the FEA and IGA results for the Tschirnhausen example under q0=250 lb/in

4.3 Tschirnhausen beam with variable curvature

The Tschirnhausen beam as a free-form curved beam with variable curvature is considered as the third example to predict the deflection results under a constant distributed load.The NURBS-based geometry of the Tschirnhausen curve is defined using the curve fitting method(see Sect.3.2).The data points for the curve fitting process of the Tschirnhausen geometry are generated using the following parametric relations:

Figure 14 illustrates the cubic B-spline curve with 30 control points representing the geometry of this curved beam.The beam,with 1 in×1 in cross section,is made of steel(E=30M psi and ν=0.25)and clamped at both ends.A uniformly distributed load q0is applied to the beam which is always in the normal direction upward.

Fig.16 Large deflection analysis of Tschirnhausen beam under uniform distributed loading:a load-deflection curve,b ratio of nonlinear to linear deflection versus load

Fig.17 Deformed geometry of the Tschirnhausen beam’s midline for q0=250 lb/in

Fig.18 Clamped-free cardioid beam under a horizontal load at the end

Fig.19 Convergence rates of the FEA and IGA results for the cardioid beam example under F=300kips(1 kips=453.6 kg)

The maximum total displacement is considered as the output of the analysis.The IGA results are obtained for two geometries with basis functions of degrees 3 and 6.This case study is also investigated in ABAQUS environment as commercial FEA software.Considering the valuable specifications of hybrid elements in incompressible materials and their superiority in avoiding volumetric locking,they are employed in the FEA solutions that smoothly improve the accuracy of results.The number of elements is well refined to reach 10?3converged outputs.Figure 15 compares the convergence rates of the FEA results of ABAQUS and IGA results with cubic and sextic basis functions for q0=250 lb/in.Since an analytical solution does not exist,the deflection error for a specific number of elements is determined as the difference between the computed displacement and the final converged value.The figure shows that the IGA approach(especially that of the sixth degree)converges faster and may need fewer elements to achieve desirable results.

Fig.20 a Load-deflection curve of the cardioid beam example;b deformed geometry of the midline of the cardioid beam under F=300kips

The comparison of results of large deflection analysis is shown in Fig.16where we can find that the presented method of this article is well consistent with the commercial FEA software.The figure shows the nonlinear load-deflection curve as well as the ratio of nonlinear to linear deflection versus load where the latter indicates how the nonlinearity increases with the applying load.Figure 17 also plots the deformed geometry of the beam’s midline for q0=250 lb/in.

4.4 Cardioid beam

The cardioid beam(Fig.18)as another free-form geometry with variable curvature is considered as the fourth case study to predict the deflection results under a point load at the end.The NURBS-based geometry of the beam is defined using the data points generated by the parametric relations of Eq.(35).The midline arc-length of the beam is 80 in and,as depicted in the figure,it represents a thick beam so that the Timoshenko beam theory should be applied.

A convergence testis also performed for the FEA and IGA outputs under F=300 kips where the elements are refined to reach the 10?3accuracy in results.As demonstrated in Fig.19,the IGA approach converges faster,especially when the 6-th degree elements are employed.The reason is that elements with the shape functions of the higher degree can better approximate the curved geometry of the beam.The nonlinear load-deflection curve and deformed shape of the midline of the beam are shown in Fig.20,assuming that the stresses remain in the elastic range.

Fig.21 Simply-supported logarithmic spiral beam under a bending moment at the top end

Fig.22 Convergence test of the FEA and IGA results for the logarithmic spiral beam under M=6 ×104 lb·in

Fig.23 a Load-deflection curve of the logarithmic spiral beam under a bending moment at the top end;b deformed geometry of the midline of the logarithmic spiral beam under M=6 × 104 lb·in

4.5 Logarithmic spiral beam

The logarithmic spiral beam is the last case study of this article that is also a free-form geometry with variable curvature.The geometry of the beam is parametrically defined by Eq.(36)in the x-y coordinate system.The curve fitting algorithm should then be employed to generate the NURBS-based geometry for IGA.As depicted in Fig.21,the simply supported logarithmic spiral beam is made of steel tube and supposed to be under a bending moment at the top end.The beam with the midline arc-length of 148 in has an overall height of 76 in and width of 43 in and thereby can be considered as a thick beam.

An important point to consider is that in the framework of isogeometric beam elements,end loads can be included in the load vector as boundary terms.The Timoshenko beam element has an interesting superiority over the Euler–Bernoulli beam element that the end moments are also emerged in the load vector since the beam has a rotational DOF as illustrated earlier in Fig.3.For the logarithmic spiral beam of this case study,the convergence test,the nonlinear load-deflection curve and the deformed shape of the beam under a constant bending moment of6×104lb in are demonstrated in Figs.22 and 23,respectively.It should be pointed out that in this case study,similar to two previous examples,the cubic elements are preferred in the isogeometric analysis despite having a faster convergence for 6-th degree elements.As a matter of fact,the degree of the NURBS curve is set to three that is the most widely accepted value in CAD applications[6].It is also worth noticing that a slight difference between the IGA and FEA results of ABAQUS might be due to the fact that commercial FEA software uses linear or quadratic shape functions for the model representation while the free-form curved geometry is of a higher degree in the NURBS format.

4.6 The parameterization study

An interesting topic in constructing free-form geometries from input data points is evaluating the effect of parameterization strategy.Generally,data points might be arbitrarily distributed and parameter assignment to these data can,therefore,influence the Jacobian and kinematics of the problem[26].As introduced in details in Sect.3.2,three frequently used parameter selection schemes are the uniformly-spaced,chord-length,and centripetal parameterization.The uniformly-spaced method does not take the initial distribution of data points into consideration,while the chord-length and centripetal methods can take care of it.Curves obtained by these three parameterization strategies are almost geometrically the same,but parametrically different and as a result,kinematic factors such as Jacobian,which are relied on the curve’s derivatives,are totally related to the parameterization[41].The parameterization is in fact a hidden concept of the IGA approach that should be brought to attention.Inasmuch as the Jacobian is frequently emerged in the IGA formulation of free-form curved beams,the influence of parameterization is investigated for three free-form examples:Tschirnhausen,cardioid and logarithmic spiral beams.The distribution of data points and deflection error with respect to the number of elements for these examples are depicted in Figs.24 and 25,respectively.

Fig.24 Data point distributions of the midline of a Tschirnhausen,b cardioid,and c logarithmic spiral beams

Fig.25 Effect of parameterization on the IGA solution of a Tschirnhausen,b cardioid,and c logarithmic spiral beams

Although all methods in all examples have been led to the same converged values,it is clear that the convergence rate in the case of chord-length procedure is higher than other two methods.In addition the effect of non-uniform data point distribution will be better reduced by chord-length parameterization[26].Hence,it is suggested to employ the chord-length method as the parameterization strategy in curve fitting.

5 Conclusions

With a wide variety of applications of curved beams with variable curvature,the objective of this article was to extend the isogeometric analysis approach to account for nonlinear force–displacement relations of free-form curved beams based on the Timoshenko beam theory.Based on the presented formulation and results,the following conclusions can be drawn:

(1)The large deformation of straight beams was successfully analyzed with IGA.The superiority of the Timoshenko theory compared to the Euler–Bernoulli theory in thick beams was investigated and showed a good correlation with other available numerical solutions.In addition,the results show that for the large applied loads under which the nonlinearity increases,the IGA solution can bring in accurate results.

(2)The validity of the model was successfully tested in the well-established shallow arch example,which represents a curved beam with constant curvature,and the Tschirnhausen,cardioid and logarithmic spiral curved beam that represent curves with variable curvature.

(3)The convergence study,which was performed for the large deformation analysis of free-form beam examples,shows that the IGA approach with basis functions of the higher order converges faster than the FEA solution since it can express a more accurate representation of the curved geometry of the beam.

(4)The parameterization study was performed to assess the parameterization effect on the IGA results for free-form geometries.The results showed that the chord-length parameter selection method converges faster and should be employed for generating suitable geometries from input data points.

As future works,it is suggested to extend the current IGA formulation of this article to spatial beams with large deformations.Moreover,it can be suggested to analyze the nonlinear dynamic behavior of Timoshenko beams with large deformation.