A FULL NONLINEAR CURVED BEAM ELEMENT WITH MIXEDFORMULATION IN STATE SPACE

In this paper, a force-based curved beam element is presented for geometrically nonlinear quasi-static analysis. The development is based on the Reissner’s exact stress resultant theory and its finite strain field for shear deformable curved beams and arches. The presented technique is found by the flexibility based method in which force interpolation functions are used. The state space approach, where a differential-algebraic equation system is solved simultaneously, is utilized as the system solution procedure. In order to improve the element accuracy, a higher order displacement field approximation is utilized based on Lagrange polynomials to evaluate the element flexibility matrix. Finally, the proposed method is validated by nonlinear examples which include some high nonlinear responses as snap backs and steep downward slopes as well as curvilinear beams and shear deformations effects. The comparison of this mixed technique with general displacement-based finite element approach demonstrates some improvement in the accuracy and reliability of the presented formulation with less discritizations. Besides, the shear/membrane locking is alleviated by the element because of using a mixed technique.


INTRODUCTION
Generally, exact and efficient analysis of frame structures, using robust numerical methods should be based on proper nonlinear beam theories such as Reissner's finite strain field [1,2] which itself is based on Timoshenko's plane cross section assumption.Since the constitutive relations are written for stress and strain resultants instead of stress and strains, the Reissner's beam theory reserves the computational efficiency in comparison with Timoshenko element model.
Many researchers have used this strain field or developed it for nonlinear analyses.Simo [3], Simo and Vu-Quoc [4] formulated and implemented an exact theory for three dimensional shear deformable beam element which reduces to Reissner's in planar case.Saje [5] used the Reissner's kinematic relations to present a finite element formulation for the finite deformation of arbitrary curved extensible shear deformable beam.Pajunen [6] performed large deflection elasto-plastic analyses implementing Reissner's kinematically exact beam theory.Jayachanrdan and White [7] implemented these strain displacement relationships for planar beam using a variable order secant matrix.
As another essential aspect of the modeling, two main points of view through the curved beam simulation have been regarded by researchers.First is using straight beam elements based on straight beam theories which generally require more elements to obtain satisfactory results., Ibrahimbegovic et al. [8]).Second is the assumption of curvilinear reference line for the curved beam based on slender beam theories to get more accurate results with lower number of elements.(Saje [5], Simo el al. [9] and Ibrahimbegovic [10]).However, the applications of these methods result in some troubles for example in long term dynamic problems including thick or moderately thick beam elements.
For the sake of completeness, it's noted that the early attempts to develop shear deformable curved finite elements were not very successful and the resulted elements might exhibit an excessive bending stiffness under inextensible deformation state and excessive shear stiffness particularly in very thin beam cases.To develop a new shear deformable curved element, it's vital to properly conquer the locking problem for both moderately thick and very thin problems.There are generally four ways to overcome this trouble [11]: (i) The assumed strain technique [12], (ii) The reduced integration scheme [13,14], (iii) The special hybrid/mixed elements [15], (iv) Appropriate Kirchoff/Mindlin representation with higher order coupled displacement rotation field [16].
Up to now, traditional computer structural analyses have generally implemented stiffness based elements that have some drawbacks associated with displacement field interpolations.This particularly happens in highly nonlinearities and for non-prismatic beams [17].Due to the time consuming solution processes in displacement based approaches, the equilibrium-based elements have been employed to enhance the procedures by strict equilibrium satisfaction along them.For comparable accuracy in global and local responses, researchers have shown the robustness of these methods even in states of the softening and using lower number of degrees of freedom.The main advantage of the flexibility methods is the low independency of the model discretization to errors in comparison with usual displacement based techniques; hence a fast convergence rate is achievable.Additionally, spread plasticity can be included in this formulation inherently.
Generally, two main aspects of element state determinations have been implemented in force based applications.In the first method which was employed primarily by Spacone et al. [18,19] and Taucer et al. [20], the residual displacement field is found in an iterative procedure on the basis of updating forces and stiffness matrices.However, in the second scheme which was initially applied by Simeonov et al. [21] and Sivaselvan and Reinhorn [22,23], the state space approach is utilized.Here, a set of Differential Algebraic Equations (DAEs) are developed consisting of rate form nonlinear material constitutive models (in the form of first order differential equations) and equilibrium equations of motion.The solution technique of the DAEs system was developed firstly by Brenan et al. [24].This technique improves iterative procedures by means of a single characteristic equation in all stages of cyclic nonlinearities.Besides, for elasto-plastic analyses a smoother transition can be obtained between all parts of nonlinear branches.
In this study the dynamical state space approach is applied to a force based shear deformable thin curved beam element as a solution procedure.The formulation is founded on the flexibility based method in which force interpolation functions and principle of virtual forces are used.In order to obtain the element flexibility matrix, the strong form equilibrium and weak form compatibility equations of the system are employed based on the Reissner's strain field.It's obvious that the equilibrium equation can be found simply for a geometrically nonlinear beam element as described by Reissner [2].However, the compatibility equations of the system in State Space approach are written in rate form.
Although the presented element has two nodes, its reference line is not straight.This is because of the existence of initial conditions at internal integration points.The results of some geometrically nonlinear bench marks demonstrate the accuracy and reliability of the element with lower discritizations.Besides, the use of the mixed techniques removes shear and membrane locking in the element.Additionally, the application of the continuous form of Arc Length constraint equation to the dynamical system leads to robust responses in high nonlinear phenomena including sharp downward slopes, snap through and snap back.

ADOPTED COORDINATE SYSTEM
Generally, a planar two node finite element curved beam has six degrees of freedom in undeformed state as shown in Fig. 1a.The element stiffness matrix is singular in this system due to presence of rigid body modes; consequently a new coordinate system is adopted so as to exclude these singularities.Fig. 1b presents this coordinate system and its degrees of freedom.This presentation has been used by Carol and Murcia [25], and Neuenhofer and Filippou [26] for straight beams.The three local degrees of freedom are (1) Where vector components are shown in Fig. 1.Since the structure undergoes large deformations with small strains, the choice of reference configuration is important for establishing the governing equations.To handle arbitrarily large rotations during the analysis, the corotational transformation method is used between two coordinate systems in this work.For more details see, e.g.[27,28].

KINEMATIC HYPOTHESIS
The implemented formulation is based on the Timoshenko beam theory assuming that plane cross sections remain plane, but not necessarily normal to the deformed longitudinal axis.Besides, changes in cross section geometry are neglected.In this regard Reissner has developed a kinematically exact stress resultant theory for curved beams and arches.The finite strains in this formulation can be specialized as In which   / 0 d ds   denotes the differentiation with respect to element undeformed axial coordinates, u and v are the displacements of the x  and y  local direc-tions, respectively and  measures the angle between beam cross section and the y  axis.
For simplicity the parameters are written in simplified form, for instance   3) and ( 5) can be rewritten as Where:  

EQUILIBRIUM AND COMPATIBILITY EQUATIONS
The strain-displacement relationship leads to the following helpful relations:     Where   ,   is the coordinate of a point which was   , x y on the initial configuration before deformation.Also / 0 s ds ds   is the differentiation of current configuration with respect to the initial state.As it'll be seen next in this article, Eqs.(11) and ( 12) will be employed for the alternative displacement interpolation technique, on the way to improve the accuracy of the responses.
If the shear and curvilinear effects are neglected, the abovementioned equations can be simplified to those of Huddleston [29], as used by Sivaselvan and Reinhorn [22,23] for the straight Euler-Bernoulli beam in the state space approach.
The geometrically nonlinear equilibrium equations for the planar shear deformable beam-column element (Fig. 1b) has been derived by Reissner [1] as: (13) In which 1 l L D   .As it's seen the relation between section force components and element end forces is shown in Eq. (13).
For compatibility equations, it's noted that a rate form is considered in this technique, which makes it different to the previous matrix-based methods by Spacone et al. [18,19].The weak form compatibility equations for Timoshenko frame can be derived using the completely analogous procedure (but with shear effects) to that of Sivaselvan and Reinhorn [22] for Euler-Bernoulli beams.From Eqs. (11) and (12) we have: After some simplifications using integration by parts, a similar result is achieved.For briefness the detailed procedure is not mentioned here, but can be traced similar to Euler-Bernoulli beam element in [23].The compatibility equation is finally derived as

SECTION CONSTITUTIVE RELATION
The section in this work is assumed to be elastic and symmetric, therefore its constitutive relation for a geometrically nonlinear element is determined by Where E T and G T are the axial and shear material tangent modulus.Also A , I and A s are section area, inertia and shear area, respectively.Eq. ( 18) has been written based on: In which       are the axial and shear stresses over the section, respectively.The symmetric section flexibility matrix is then computed as the inverse of section stiffness matrix as follows:

FLEXIBILITY-BASED PLANAR CURVED BEAM ELEMENT FORMULATION
The flexibility matrix for a geometrically nonlinear is obtained using the rate form compatibility equation.In this regard Eq. ( 17) can be rewritten as Therefore the element flexibility matrix is It's noted that the first term (at the second side) in Eq. ( 21) indicates the initial stress term.Both Gauss and Gauss-Lobatto rules are valid for the flexibility matrix evaluation, but as a more applicable method for general structural analysis, the Gauss-Lobatto rule is utilized to evaluate the flexibility matrix in the abovementioned equations.

APPROXIMATION OF THE ELEMENT DISPLACEMENT FIELD
While the element flexibility matrix is evaluated at Gauss points inside the element, the element displacement fields including  ,  and  need to be evaluated at these points.Additionally, the element accuracy in high nonlinearities can be directly affected by the approximation procedure.That's why a higher order displacement approximation method (in comparison with that of Sivaselvan and Reinhorn [22]) is adjusted in this work.To this end a polynomial based approximation field is developed for the beam element as described hereafter.
The proposed technique is basically a combination of linear approximation field [23] and curvature based displacement interpolation (CBDI) technique [26].In the present formulation, the deformations are monitored at NG integration points in terms of NG parameters as   , where j s are the local Gauss points coordinates on the initial curved beam configuration.According to this: Where ( ) l s j are the Lagrange Polynomials which are obtained from the following formula While there is a need to integrate over functions, the more applicable form of the polynomials for closed form integration is used as Where G is the so called Vandermode matrix described by As the number of integration points doesn't change during the analysis, this matrix is determined only once.Starting from one end of the element, the incremental form for the internal displacement field can be stated as To use a higher order and more accurate integration, the incremental values are found as Where 0 S indicates the undeformed length.Once, the displacement field is obtained from one element end, the procedure is repeated from the other end and the same operations as Eqs.( 30) through (35) are performed.Finally a weighted average of two resulted fields is taken and three required parameters are obtained.In this work, the weights have been considered equal for every internal node.Besides, the local boundary conditions are enforced to end nodes.
The integrations in Eqs. ( 34) and ( 35) are achieved based on a local Gauss or Gauss-Lobatto procedure which is performed between s s i  and s s j  .

SUMMARY OF STATE EQUATIONS USING CONTINUOUS SPHERICAL ARC LENGTH TECHNIQUE
For the sake of completeness, it's noted that a general form which has been proposed before [22], is taken into account for the global and local state equation sets.However, in order to capture the snap back or steep negative slopes in responses during a quasi-static analysis, a spherical Arc Length technique is used in this study to develop the continuous constraint equation.In this state, the constraint equation is established through the controlling response arc length rather than single characteristic displacement controlling methods or load incremental controlling procedures.For more details see e.g.[30,31].The simple form for the additional constraint equation is where s is the incremental arc length.Also   f t is a function which demonstrates the direction of moving in the quasi static analysis.u  and P  are the derivation of displace- ment and applied force with respect to t , respectively.Besides,  is a scaling parameter for loading terms.Consequently, the whole form of the DAEs set is described as follows: -Nodal equilibrium equations (Three equations at each node for a two-dimensional problem) -Boundary conditions (specified displacement components at some boundary nodes) -Section constitutive relations (Including three sets of differential equations for axial strain, curvature and shear strain in every element) -Element state equations ( six equations exhibiting the relationship between end forces and end displacement components for each element)

Unrolling And Re-Rolling of A Semi Arch Beam With Moment At The Tip
The first example represents the geometrically nonlinear analysis of a semi circular cantilever beam with a bending moment applied at the free end.This demonstrates not only the ability of the co-rotational formulation to solve finite strain problems (with a suitable media discritizations) but also the efficiency of the method to handle large initial curvatures in the beam.This problem which its configuration and properties are shown on Fig. 2, has been solved previously by Jayachandran and White [7].The beam has a rectangular cross section of B H  .The quasi-static analysis is performed using five elements and five internal Lobatto points.As shown in Fig. 2 the beam deforms form the initial semi arch state to direct, semicircular arch and then a full circle in the reverse situation.Table I shows the results of the analysis for vertical and horizontal displacements of the tip using the proposed element for different moment values.It's observed that using five elements in this problem results in approximately 1.65% relative error in end rotations.

Figure.2 Large deflections of a flexible bar subjected to end moment It's noted that the use of Lobatto points instead of Gauss points is because of the more applicability of these type in the elasto-plastic structural analyses i.e. existence of two end nodes in Lobatto scheme.
Also, the values of shear forces and shear deformations in this problem are found almost exactly as zero.It's seen that while these values are independent from the height of the beam element, no shear/membrane locking has been occurred during large deformations.This is because of the utilization of a force based method, which itself is a particular state of mixed techniques.The next example provides more detailed information on the removing shear/membrane locking effect.

A Tip-Loaded Cantiliver Ciecular Arch
This example reveals that the shear/membrane locking phenomena is alleviated in the present formulation.A quarter ring subjected to radial point load P at its free end is shown in Fig. 3.The exact solutions for three displacement vector components at the tip are mentioned using the Castigliano's energy theorem in [11].To study the shear and membrane locking effect in the presented formulation, a range of / R h ratios between 50 and 1000 are utilized and the results are mentioned in table II.It's seen that both shear locking and membrane locking are removed for different values of / R h ratios for very slender states and the results are almost independent from the / R h ratios.Consequently, the formulation is reliable to be applied to very thin beams without any problem.

Figure.3 A quarter circular cantilever ring
Table II.Displacement components of the tip of quarter cantilever beam subjected to end load normalized to exact solutions [11].u : horizontal displacement, v : vertical displacement and  : end rotation.

Post-Buckling Analysis of A Clamped-Hinged Circular Arch
This well-known example presents the geometric nonlinear analysis of the pre-and post-buckling deformation of a circular arch, hinged at one end and clamped at the other end, under a vertical force applied at the apex.This problem was suggested by Da Deppo and Schmidt [32] and has been analyzed by many researchers.(Simo and Vu-Quoc [4], Jayachandran and White [7] and Wood and Zienkiewicz [33] , Ibrahimbegovic [34]).
Material and geometric data for the arch are shown in Fig. 4. As responses show, the problem includes a sharp negative slope after the first buckling load is reached.The use of continuous form of constraint equation results in robust responses as are demonstrated in Fig. 5 for eight elements and five Lobatto points.Also, table III compares the results of the proposed technique (using eight and ten elements, both with five internal Lobatto points) with previous researches.

Model
Num. of Elems.

Post-Buckling Analysis of A Shallow Circular Arch
A shallow circular arch which is pinned at its two ends and is subjected to a concentrated load at the point close to its peak is shown in Fig. 6.The properties of the problem have been illustrated on the figure.This bench mark problem has been analyzed by some researchers such as Li [35] and Clarke and Hancock [36].For this problem two elements with five Lobatto points inside each element are considered.The load deflection response curve of this shallow arch is depicted in Fig. 7, where multiple load critical points and also displacement limit points exist at the buckling and post buckling stages.The fitness of the present study with the results in several critical points obtained by Clarke and Hancock [36] and Li [35] has been demonstrated in Fig. 7. It's seen that since the exact equilibrium is satisfied along the curved elements, accurate and reliable results have been found with the present flexibility based technique using less discritizations.The results of the present method are comparable with those of Li [35], with 11 elements.Some deformed states of the arch are shown in Fig. 8.

CONCLUSIONS
A formulation is presented for a geometrically nonlinear force-based shear deformable curved beam element.The curved inherent of the presented two node element is hidden in its integration points.In this approach the elements geometrically nonlinear behavior is formulated by differential equations involving a set of variables that describe the state of each element and solve them simultaneously with the global equations of motion.While the presented formulation is based on the mixed approach, it removes shear/membrane locking inherently.Additionally, the application of the new displacement approximation filed with the force based technique in which equilibrium equations are satisfied exactly along the element, leads the analysis to more accurate and reliable results with less discritizations.Besides, the use of this presentation in conjunction with the continuous form of the Arc Length additional constraint equation makes it possible to trace the complete equilibrium path which may include sharp downward slopes or snap backs, very stably.The presented technique is very straight forward to be applied to elasto-plastic analyses, because a single characteristic equation in all stages of cyclic nonlinearities is used.Moreover, the effects of interaction can be incorporated in the analyses, using the same procedure described by Sivaselvan and Reinhorn [22].

Figure. 1 .
Figure.1.Global and Local coordinate systems, (a) Global coordinate system and degrees of freedom for a planar undeformed curved beam element.(b) Local adopted coordinate system without rigid body modes in element deformed state

Figure. 4 .
Figure.4.Circular arc geometry and deformed shape in some different load levels

Figure. 6 .
Figure.6.shallow circular arch subject to an offset concentrated load

Fig. 8
Fig.8 the deformed shapes of shallow circular arch under different load levels

Table I .
Horizontal and Vertical displacements of the tip of semi circular cantilever beam subjected to end moment compared with analytical solution which is found regarding,