Developing a Formulation Based upon Curvature for Analysis of Nonprismatic Curved Beams

A new element with three nodal curvatures has been considered for analysis of the nonprismatic curved beams by finite element method. In the formulation developed, the force-curvature relationships in polar coordinate system have been obtained first, then the curvature of the element has been assumed to have a second-order polynomial function form and the radial, tangential displacements, and rotation of the cross section have been found as a function of the curvature accounting for the effects of the cross section variation. Moreover, the relationship between nodal curvatures and nodal deformations has been calculated and used for determining the deformations in terms of curvature at an arbitrary point. The total potential energy has been calculated accounting for bending, shear, and tangential deformations. Invoking the stationary condition of the system, the force-deformation relationship for the element has been obtained. Using this relationship, the stiffness matrix and the equivalent fixed loads applying at the nodes have been computed. The results obtained have been compared with the results of some other references through several numerical examples. The comparison indicates that the present formulation has enough accuracy in analysis of thin and thick nonprismatic curved beams.


Introduction
Because of much application of curved beams and arches in different structures, in recent decades many studies have been conducted for analysis of such beams using finite element methods, Marquis and Wang [1], Litewka and Rakowski [2,3], Friedman and Kosmatka [4].Finite element analysis of the curved beams using several straight beam elements, or low-order curved beam elements have been proposed by Mac Neal and Harder [5]; but the effects of shear and membrane locking exist in their results.Ashwell et al. [6,7] have analyzed the curved beams using two models; in the first model the radial deformation has been approximated by a third-order polynomial function, but the tangential deformation by a linear function, and in the second model the shape functions for the cylindrical shell elements have been used in order to determine the deformations of the curved beam elements.Comparing the results of two models has indicated that the second model has been more accurate for thick arches with low displacements.Dawe [8] has introduced eight different models using different shell element theories.Among the models introduced, the results obtained by using high-order shape functions are more accurate even though modeling them is more difficult.Krishnan and Suresh [9] have calculated the deformation of an arch element by static and dynamic analysis of the arch using a pair of polynomials of order three that have been guessed, and then they have discussed the error resulting from eliminating the effect of the shear in determining the natural frequency of the arches.Babu and Prathap [10] have introduced the radial and tangential displacements and rotation for several models of the thick curved element guessing linear shape functions.In these models, the error arising from the shear locking has been emphasized.The relationships for the curved beam element have been determined without the effect of shear using penalty method by Tessler and Spiridigliozzi [11]; but the relationships of the proposed model are so complex.Stolarski and Belytschko [12] have considered the membrane locking phenomenon in thin arches with low displacements using high-order polynomials.The relationships obtained have a significant error in the analysis of thick arches because of shear and membrane locking.
In recent years, the locking phenomenon has been paid much attention by researchers.When the element is under bending and elongation of the fibers is restricted, the membrane locking will be produced.On the other hand, if in the element formulation, the effects of pure bending are considered without implementing the effects of shear, then the shear locking will be existed.Lee and Sin [13] have presented a three-node element based upon curvature having constant cross section.They have studied the results of the deformations for different thickness values through numerical examples.The results obtained indicate that the shear and membrane locking phenomenon has been reduced significantly.Raveendranath et al. [14] have introduced a two-node arch element accounting for shear effects, but they have reported that the results for high thicknesses are with error.In several references ignoring the membrane effects, the effects of shear in static analysis of prismatic curved beams have been discussed (Yang and Sin [15], Sheikh [16]).Chen [17] has studied the effects of shear for loading inside the curved beam plane using cubic technique.In the previous studies carried out by Sinaie and the authors of this paper [18] considering the relationships proposed by Lee and Sin [13], the high-order curvature shape functions have been used instead of displacement shape functions and then, the force-deformation relationships have been determined.Sinaie et al. [18] have eliminated the shear and membrane locking in the relationships proposed by Lee and Sin [13] for prismatic curved beam.In other words, in their work the effects of the shear and tangential deformation and also the effect of the bending moment in the calculation of the curved element with six nodal curvatures for analysis of thin and thick prismatic beams have been accounted for, which is the main difference with the method proposed in other references.In the present paper, after finding the curvature-strain relationships and the internal forces in polar coordinate system, guessing a second-order polynomial function for the element curvature, the shape functions have been determined first.Then, the curvature at any point of a three-node element has been found based upon the nodes curvatures.Furthermore, the relationship between the nodal curvature and the nodal displacement has been found accounting for the effects of the cross section variation using a transform matrix.The ratio of the moment of inertia to the cross sectional area has been considered as a function of the arch length.Finally, minimizing the total potential energy, the forcedeformation relationship and the stiffness matrix in local coordinate system have been found.Moreover, an algorithm for analyzing the nonprismatic curved beams has been presented.At the end, through three numerical examples, the results of using the method proposed in this work have been compared with the results of (a) using shell elements (b) using prismatic curved beam elements (c) exact solution.It has been indicated that since in the proposed element the shear and membrane deformation effects have been accounted for, the accuracy for analysis of thin and thick nonprismatic curved beams has been increased.The simplicity of modeling and high accuracy of the analysis are the privileges of the proposed method compared to the methods of the other references.In the other methods, the member must be divided into many elements, while in the proposed approach the convergence will be achieved with one or two elements only.

Formulation of the curved element based upon curvature
The formulation for an element of the nonprismatic curved beam shown in Figure 2.1 will be conducted in four steps.The first step is setting up the basic curvature-displacement equations for an element of the curved beam in polar coordinate system and determining the components of the radial, tangential displacements, and rotation in terms of curvature considering the variation of the cross sectional area; the second step is determining the shape functions that state the curvature of any point on the element axis in terms of the curvatures of three nodes adopted, which is done by choosing a secondorder function; the third step is finding a relationship between the nodal curvature at the three nodes and the nodal displacements at two ends of the element; the fourth step is calculating the total potential energy considering the effects of the internal membrane and shear forces and bending moment.Invoking the stationary condition of the system, the force-displacement relationship, and then, the stiffness matrix of the member and the equivalent fixed load vector can be found.

2.1.
The basic curvature equations in polar coordinate system.In Figure 2.1, an element of nonprismatic circular curved beam having radius R and the arch length L is shown.The displacement components are the radial displacement, W, the tangential displacement, U, and the rotation, θ.The strain-deformation relationships are as follows (Timoshenko and Goodier [19]): ) ) where ε is the tangential strain, κ is the curvature, γ is the shear strain and the subscript ,s indicates the derivative with respect to the arch length.
The relationship between the internal forces of the element and the strains can be stated as follows where N is the tangential force, M b is the bending moment, V is the shear force, A and I are the cross sectional area and the moment of inertia as functions of arch length s, respectively, n is the shear correction coefficient, E and G are elasticity and shear modulus, respectively, and κ is the curvature of the deformed beam.On the other hand, setting up the equilibrium at point s of the element, the following relationships are obtained: Substituting (2.4) into (2.5) and (2.6), the relationships between strains and curvature can be found as follows: (2.8) H. Saffari et al. 5 Substituting (2.8) into (2.1) to (2.3), the displacement-curvature relationships are obtained as follows: ) From (2.2), the rotation θ in terms of the curvature κ can be calculated as follows: in which the constant C θ can be computed based upon the end conditions.Moreover, the following differential equation which states W in terms of κ will be obtained by eliminating U in (2.9) and (2.10): (2.12) The solution of the differential equation (2.12) consists of general solution, W h , and particular solution, W P , as follows: ) The constants C W1 and C W2 will be calculated applying the end boundary conditions.W P depends on A, I, and κ which all are functions of arch length s.Substituting W from (2.13) into (2.9), the tangential displacement U can also be computed as follows:

Determining the relationship between the curvature at an arbitrary point and
Nodal curvatures.In this section, guessing a second-order polynomial function for curvature and conducting some regular matrix operations, the shape function and also the curvature function, κ, at any point will be obtained.Therefore, the relationship of the element curvature at any arbitrary point in terms of the curvatures of three nodes will be found (see Figure 2.2).As indicated in Figure 2.2, the curvature has been considered at three points of the element, so, the curvature function has three unknown constants.Hence, Matrix form of the function (2.16) arises from multiplying geometrical vector [g] by constant coefficients vector [C] as follows: (2.17) If Q is the vector indicating the curvature at three nodes as shown in Figure 2.2, then where g i is g at point i.So, the relationship between the nodal curvature and the curvature at an arbitrary point is stated as follows: ) The components of the shape functions vector concerning the curvature function, F κ , are shown as follows: The displacement components of the element (radial, tangential, and rotation) will be computed by substituting the curvature function from (2.20) into (2.11),(2.13), and (2.15) as follows, ) (2.28) F κ , F θ , F W , and F U can be found having A and I functions.The values of these functions for the case that the cross section of the beam is rectangular such that its height is a function of arch length s can be found in Appendix A.

Determining the relationships between the Nodal curvatures and the Nodal displacements.
Using (2.23), (2.25), and (2.27), the components of the radial and tangential displacements and rotation at node 1 of the curved beam element with 6 degrees of freedom, as shown in Figure 2.3, are (2.31) The components of the radial and tangential displacements and rotation at node 2 of the element are (see Figure 2.3) ) ) If (2.35) are rearranged in the matrix form, the following relationships between nodal displacements and nodal curvatures will be obtained: where (2.37) The value of the curvature at any arbitrary point of the element in terms of nodal displacements can be calculated as follows: displacements can be computed as follow: ) ) in which F W | 0 and F U | 0 indicate the values of F W and F U at s = 0, respectively.The values of F θ0 , F W0 , and F U0 can be found in Appendix B.

The equilibrium equation of the finite element and force-deformation relationship.
The total potential energy of a curved element as shown in Figure 2.2 accounting for the effects of tangential and shear forces and bending moment can be presented as follows: in which r and t are the radial and tangential distributed external loads, respectively, and m is the distributed bending moment as indicated in Figure 2.2.Furthermore, U e is the strain energy which can be found as follows: The components of the radial and tangential displacements and the rotation from (2.39), (2.40), and (2.41), respectively, will be substituted into (2.8)first, and then, the result will be inserted in (2.43).So, accounting for the effects of the membrane, shear and bending deformations in the total potential energy, the shear and membrane locking will be eliminated and it is expected that such element can be used for analysis of the thick beams besides the thin beams.Now, the stationary condition of the system, δΠ = 0, is invoked using Mathematica Software [20].and the result is arranged in the force-deformation relationship as follows: in which [K] is the stiffness matrix and can be found as follows: where f (s) = I/A (the ratio of the moment of inertia to the cross sectional area as a function of the arch length, s).
[P] is the resultant of the forces acting on each node which can 10 Mathematical Problems in Engineering be found from the following relationship: The values of S W , S U , and S θ can be found in Appendix B.

Numerical studies
In this section, the formulation presented in this paper has been applied to three nonprismatic curved beams with different characteristics and different end conditions, and the results obtained have been compared to the results of the finite element model with shell elements, to the results of the energy method, and also to the results presented by Lee and Sin [13].The first example concerns a curved beam in the form of a quarter of a ring having one end clamped and a concentrated load applied on the other end.

Example 1: a quarter of a ring.
A curved beam in the form of a quarter of a ring, as shown in Figure 3.1, has been considered for Example 1.Its fixed end is setting on the origin of the coordinate system and the external load is applying on its free end.In this example, the cross section of the beam is rectangular having width equal to unity.The H. Saffari et al. 11 height of the cross section, H(s), is a linear function of the arch length, s, as follows: ) where H1 is the height of the section at its fixed end and L is the length of the curved beam.This example has been solved for different ratios of H1/R in order to indicate the capability of the formulation developed which considered the high effects of shear deformation for thin and thick beams.For better comparison, four methods have been used for solving this example as follows.
In the first method, the displacements have been found using Castigliano's theorems through the following processes.
The internal forces of the ring are On the other hand, the strain energy of the curved beam element is 12 Mathematical Problems in Engineering Radial displacement can be obtained using Mathematica Software [20] by taking ds = R • dθ and derivation of the strain energy as Then, the radial displacement based upon Castigliano's method without accounting for shear and membrane effects have been calculated as follows: The redial displacements by Castigliano's method accounting for shear and membrane is divided by the results obtained from Castigliano's method which does not account for shear and membrane effects.These ratios are indicated by W C 2 / W C 2 and shown in Table 3.1.
Similarly, the tangential displacement and the rotation have been computed and the results have been divided by the results obtained using Castigliano's method which does not account for shear and membrane deformations effects.These ratios are indicated by U C 2 / U C 2 for tangential displacement and by θ C 2 / θ C 2 for rotation.In the second method, the ring has been modeled by one element based upon the formulation developed in this paper and then, the displacements obtained have been divided by the amounts obtained using Castigliano's method which does not account for shear and membrane deformations effects.These ratios are indicated by W F 2 / W C 2 for radial displacement, by U F 2 / U C 2 for tangential displacement, and by θ F 2 / θ C 2 for rotation.In the third method, the ring has been modeled and analyzed by ANSYS [21] using one hundred elements.The number of elements has been taken high in order to increase the accuracy for the thick elements.Then, the displacements obtained have been divided by the amounts obtained using Castigliano's method which does not account for shear and membrane deformations effects.These ratios are indicated by W S 2 / W C 2 for radial displacement, by U S 2 / U C 2 for tangential displacement, and by θ S 2 / θ C 2 for rotation.In the forth method, the ring has been modeled and analyzed using 20 prismatic curved beam elements as proposed by Lee and Sin [13].Then, the displacements obtained have been divided by the amounts resulting from Castigliano's method which does not account for shear and membrane deformations effects.These ratios are indicated by W CB 2 / W C 2 for radial displacement, by U CB 2 / U C 2 for tangential displacement, and by θ CB 2 / θ C 2 for rotation.As explained, for better comparing and for considering the shear locking and membrane effects, the amounts obtained using all the four above-mentioned methods have been made dimensionless dividing by the amounts obtained using Castigliano's method H. Saffari et al. 13  The superscripts C, F, S, and CB indicate the amounts obtained using Castigliano's method, the method developed in this work, using shell element of ANSYS [21] and using the element developed by Lee and Sin [13], respectively.
θ C 2 are the amounts of radial displacement, tangential displacement, and rotation resulting from Castigliano's method which does not account for shear and membrane effects, respectively.

H1/R
Radial displacement ratio of node 1 Radial displacement ratio of node 2 The superscripts C, F, S, and CB indicate the amounts obtained using Castigliano's method, the method developed in this work, using shell element of ANSYS [21] and using the element developed by Lee and Sin [13], respectively.
W C 1 , W C 2 are the amounts of radial displacements of nodes 1, 2 resulting from Castigliano's method which does not account for shear and membrane effects, respectively.which does not account for shear and membrane deformations effects, and have been shown in Table 3.1.So, the effects of the shear and membrane phenomenon on the results can be easily evaluated and compared.
As it can be seen in Table 3.1, in thin beams, because of shear and membrane deformations effects, the results are close to 1, but in thick beams (i.e., the thickness to radius ratio more than 1.2) the results are so far from 1. On the other hand, the results obtained applying the method developed in this paper are so close to the results of analyzing by ANSYS [21].That is because of high effects of shear phenomenon, which is one of the abilities of the element introduced in this work for analyzing the thick beams.The other point is that in the present method which uses curved beam element with nonprismatic cross section, just one element has been used for convergence.It is assumed that the height of the section varies linearly based upon (3.1).This example has been considered for different ratios of H1/R using four methods, as explained in the following, and the results are indicated in Table 3.2.

Example 2: a pinched ring.
In the first method, the radial displacements at nodes 1 and 2 have been calculated using Castigliano's theorems.In the second method, the problem has been modeled using the formulation developed in this paper and the radial displacements at nodes 1 and 2 have been determined under concentrated load accounting for effects of variation of the height of the section based upon (3.1).In the third method, the analysis has been done by ANSYS [21] using 100 SHELL63 elements.Finally, in the fourth method, the ring has been modeled and analyzed using 20 prismatic curved beam elements as proposed by Lee and Sin [13].
The results of the four methods have been made dimensionless dividing by the displacements obtained using Castigliano's method which does not account for shear and membrane deformations effects, and compared in Table 3.2.As it can be seen, the difference between the results and unity indicates the high effects of shear and membrane deformations.The results of applying the method developed in this work have high accuracy although just one element has been used for analysis.Figure 3.3 shows the variation of the radial displacement component, W, versus the variation of the angle, φ, (see Figure 3.2), for a curved beam with thickness to radius ratio equal to 1/100 under unit load using the Castigliano's method and the method developed in this work.In Figure 3.3, the accuracy of the method developed in this paper is shown.The maximum error is less than 0.1%.

Example 3: simply supported thick curved beam.
In this example a simply supported thick curved beam is loaded as shown in Figure 3.4.The height of the section varies linearly based upon (3.1).The radial displacement at node 2 accounting for section variation under concentrated loading has been computed first using the method developed in this paper and then, the results have been compared with the results of the analysis method using SHELL63 element as shown in Table 3.3.As it can be seen, the results are very close.In this example, the support at node 2 is oblique which can be easily modeled by the method developed in this paper.

Summary and conclusion
In this paper, for analyzing the curved beams with variable moment of inertia, a new element has been proposed.In the proposed element, curvature shape functions have been used instead of displacement shape functions.The displacements of the ends of the element have been related to the beam curvature, and the variation of the cross sectional area has been used in the formulation of the element as a function of the arch length.Then, the force-displacement relationship has been obtained by invoking the stationary condition of the system.So, calculating the stiffness matrix and equivalent fixed end forces, a good approximation for linear analysis of nonprismatic curved beams has been presented.In this work, accounting for total bending, shear and membrane potential energy, the effect of the shear locking has been eliminated.As indicated by numerical examples, the proposed method can be used for analyzing both thick and thin beams.Simplicity in modeling the nonprismatic curved beams is one of the privileges of the proposed method.In the developed formulation, for the solution convergence one or two elements are needed only.The accuracy of the method is high while the degrees of freedom are low.
The second example considers a pinched ring under concentrated symmetric loading.The third example is a simply supported curved beam under a concentrated load.The three examples are shown in Figures 3.1, 3.2, and 3.4, respectively.

Figure 3 .
Figure 3.2.A pinched nonprismatic ring with a load of 2P and modeled for 1/4 circular arch.

Figure 3 .
2(a) indicates a pinched ring with nonprismatic cross section under two concentrated symmetric loads.Because of symmetry, the problem has been modeled as a quarter of a ring with the boundary conditions shown in Figure 3.2(b).

Table 3 .
1.Comparing the displacements at the free end (node 2) of a quarter of a ring (Example 1).

Table 3 .
[21]omparing the displacement of node 2 in Example 3.Superscripts F and S indicate the results using 2 curved beam elements (present work) and analysis by ANSYS[21]using shell elements, respectively.