ISOGEOMETRIC ANALYSIS OF PLANE–CURVED BEAMS

ANTONIO CAZZANI, MARCELLO MALAGÙ, AND EMILIO TURCO

Abstract. A curved beam element based on Timoshenko model and NURBS

(Non-Uniform Rational B-Splines) interpolation both for geometry and dis-
placements is presented. Such element can be used to suitably analyze plane
curved beams and arches. Some numerical results will explore the effectiveness
and accuracy of this novel method by comparing its performances with those
of some accurate finite elements proposed in the technical literature and also
with analytical solutions: for those cases where such closed form solutions were
not available in the literature, they have been computed by exact integration
of the governing differential equations. It is shown that the presented element
is almost insensitive to both membrane- and shear-locking, and that such phe-
nomena can be easily controlled by properly choosing the number of elements
or the NURBS degree.

1. Introduction
Finite element models usually do not replicate the exact geometry of a solid de-
fined by a Computer Aided Design (CAD) representation: the finite element mesh is
only an approximation of the exact CAD geometry, so the precision of the finite ele-
ment (FE) analysis can be affected by errors originating from geometric description.
Isogeometric analysis, whose name was coined by Hughes [1], attempts to avoid this
problem, by integrating the geometry model defined by CAD into the FE analysis
in order to solve the underlying problem without any shape approximation.
In the last years, many efforts in the framework of isogeometric analysis have
been done (see e.g. [2] and [3]) and there are interesting results for a wide range of
problems such as vibrations and wave propagations, nearly incompressible solids,
fluids, fluid-structure interaction. The geometric settings, however, mostly concern
2-D and 3-D continuum models or shell-like structures [4].
It appears that, up to the last two years, special structural elements for 1-D
problems based on the isogeometric concept were still lacking, even if the use of
B-spline had been reported, for instance, in [5], [6] and [7]. Recently, however, a
strong interest for such elements has been developing, and several contributions
appeared: see, for instance [8], [9], [10], [11], [12], [13], [14], [15]. In these papers
both straight and spatial Timoshenko beams have been considered, with particular
emphasis on locking control.
Here, a plane curved Timoshenko beam element based on NURBS interpolation
for both geometry and displacements, which is able to analyze plane arches and
beams, has been developed. In Section 2 we discuss the basic key points of the
mechanical model besides to the guidelines of the NURBS interpolation giving

Date: Published on Mathematics and Mechanics of Solids (2014), DOI:

10.1177/1081286514531265.
Key words and phrases. NURBS, Isogeometric analysis, curved Timoshenko beam.
1
2 ANTONIO CAZZANI, MARCELLO MALAGÙ, AND EMILIO TURCO

all the details useful to implement the numerical model. Successively, in Section 3,
various numerical results are presented with the goal to explore the effectiveness and
accuracy of this novel method by comparing its performances with those of finite
elements already available in the technical literature and with analytical solutions.
In particular, some exact solutions, which were not yet available in the literature,
have been computed in order to test the numerical solutions and verify the reference
values used by other authors for comparing their results. Finally, in Section 4, some
concluding remarks and future developments are briefly reported.

2. NURBS representation of plane curved Timoshenko beams

We consider a curved plane beam (see Figure 1) whose centroid line is a plane
curved parametrized by the arc-length s ∈ [0;
]. We suppose also that one of
the principal inertia axis and the shear center of the cross-section lie in the same
plane. Let the global reference system be denoted by (O; x1 , x2 ) and the local one
by (o; t, r), t and r being respectively the tangent and the normal unit vectors
to the curve. Moreover, we denote by R the curvature radius and by (·)
the
derivatives with respect to the arc-length s. With this notation, the differential

ϕ t, u
o
r, w

M M
x2 N N
T T
x1
O

Figure 1. Plane curved beam: centroid line and reference systems.

form of equilibrium and kinematic compatibility equations and the constitutive law
describing the curved Timoshenko plane beam problem in the local reference system
are:
T N
N
− + qt = 0, T
+ + qr = 0, M
− T + m = 0, (1)
R R
w u
ε = u
− , γ = w
+ + ϕ, χ = ϕ
, (2)
R R
N = EAε, T = GAT γ, M = EIχ. (3)
Here, N , T and M denote the generalized stresses (axial and shear force and bending
moment, see Figure 1 for the definition of positive quantities) and qt , qr and m the
generalized external forces per unit length (tangent and radial forces and distributed
couple moments). In the local reference system u and w are the displacements of
the axis line and ϕ the section rotation while ε, γ and χ denote the generalized
strains (axial, shear and curvature bending). Finally, symbols E, G, A, AT and I
indicate, respectively, the Young’s modulus, the shear modulus, the cross-section
area, the shear reduced cross-section and the area moment of inertia. The curved
 ISOGEOMETRIC ANALYSIS OF PLANE–CURVED BEAMS 3

beam problem can be set in an equivalent variational formulation, such as the

classical principle of total potential energy, which is certainly more suitable for a
solution strategy based on finite elements:

   
1  2 2 2
arg min (EAε + EIχ + GAT γ )ds − (qt u + qr w + mϕ)ds . (4)
u,w,ϕ 2 0 0

As a first step in order to discretize (4), we describe the geometry of a curved

beam by means of NURBS interpolation (for a complete and more accurate de-
scription of this kind of interpolation we refer to specialized books, [16]).
We say that x has a p-degree NURBS representation when there exist n ∈ N,
control points Pi ∈ R2 , weights gi ∈ R, i = 1...n, and a knot vector, i.e. a set
Ξ = {0 = ξ1 ≤ ξ2 ≤ ... ≤ ξn+p+1 = 1} such that, for any ξ ∈ [0; 1]:
n

x(ξ) = Ri,p (ξ)Pi , (5)
i=1

where the NURBS basis {Ri,p (ξ)} can be expressed as:

Bi,p (ξ)gi
Ri,p (ξ) = n , (6)
i=1 Bi,p (ξ)gi

in terms of B-splines bases {Bi,p (ξ)} defined by the Cox-De Boor recursive formula:

1 if ξi ≤ ξ < ξi+1
Bi,0 (ξ) = , (7)
0 otherwise
ξ − ξi ξi+p+1 − ξ
Bi,p (ξ) = Bi,p−1 (ξ) + Bi+1,p−1 (ξ). (8)
ξi+p − ξi ξi+p+1 − ξi+1
The so-called knot vector Ξ defines a partition of the parameter space [0; 1] simi-
lar to the classic finite element subdivision. Non-uniform knot vectors and repeated
knots are the key ingredients of NURBS flexibility and produce refined geometric
descriptions. Weights gi related to i-th control point enlarge the capabilities of the
B-splines interpolation allowing also an exact representation of conic sections.
Among all the properties of NURBS interpolation the most interesting is the
high-degree of continuity. More precisely, each p-th order function is of class C p−1 ,
i.e. it is continuous with its derivatives up to the (p − 1)-th order, and in particular
it is smooth across the boundaries. However, if necessary, continuity degree can
also be lowered by using repeated knots.
The main idea of the isogeometric approach is to exactly describe the geometry
of the problem by NURBS interpolation and to use the same interpolating basis to
represent the generalized displacements:
n
 n
 n

u(ξ) ≈ Ri,p (ξ)ui , w(ξ) ≈ Ri,p (ξ)wi , ϕ(ξ) ≈ Ri,p (ξ)ϕi , (9)
i=1 i=1 i=1

by means of control points ui , wi , and ϕi .

Using the dot to denote the derivatives with respect to ξ and denoting by J the
Jacobian of the transformation, we have:

J3
J = ṡ = ẋ21 + ẋ22 , R= , (10)
|ẋ1 ẍ2 − ẍ1 ẋ2 |
4 ANTONIO CAZZANI, MARCELLO MALAGÙ, AND EMILIO TURCO

thus, the generalized strains take the form:

u̇ w ẇ u ϕ̇
ε= − , γ= + + ϕ, χ= . (11)
J R J R J
By using the position:
 ξe+1
1
πe = (EAε2 + EIχ2 + GAT γ 2 )Jdξ,
2 ξe
 ξe+1
(12)
λe = (qt u + qr w + mϕ)Jdξ,
ξe

problem (4) can now be discretized as:

n 
e

arg min πe − λe . (13)

u,w,ϕ
e=1

where ne is the number of subdivision of the parameter space. Imposing (13)

produces a linear system of equations where the unknowns are the control points
ui , wi e ϕi .
To set up the stiffness matrix and the load vector it is necessary to evaluate
the integrals: this can be performed numerically by using the Gauss quadrature
rule, even if the integrand functions are somewhat different from polynomials. By
referring to [3] for some guidelines about the number of Gauss points to be used for
an efficient quadrature, here we resort to a simple rule suggested by some numerical
tests and validated by the numerical results presented in the next section. This is
also useful to prevent locking problems such as described in [15].

3. Numerical results
In this section we present a series of numerical tests devoted to explore the
performances of the proposed isogeometric finite element family. Their results are
compared with analytical solution, classical Lagrangian compatible elements, and
some special (i.e. mixed or hybrid) finite elements reported in technical literature.
Comparison with reference solutions is particularly devoted to:
• asses the convergence rate of the proposed element also by comparison with
classical elements such as those belonging to the Lagrangian family;
• test the influence of locking phenomena (if any);
• check these results against those produced by special elements particularly
effective for curved beams;
• explore influence of the Gauss point number on the numerical results;
• make a critical analysis of the results produced by the NURBS element
family devoted to select the most effective order.
All the numerical tests are built by using an ad hoc developed code which is
essentially based on NURBS algorithms reported in [16] and coded in a NURBS
library and on the GeoPDEs library of [17].

3.1. Tests devoted to calibrate the numerical model. Following the guide-
lines developed in [3] and bearing in mind the observation reported in [15], we
consider the simple numerical tests reported in Figure 2. For both tests a reference
solution is easy to evaluate. Furthermore, the first is a pure bending problem and
 ISOGEOMETRIC ANALYSIS OF PLANE–CURVED BEAMS 5

the second a pure extensional problem. So they can use to highlight possible locking
problems.

q
R

(a) Pure bending (b) Pure extension

Figure 2. Sketch of cantilever circular arch subjected to a couple

at the free end and of ring under internal pressure.

First, we consider the cantilever circular arch sketched in Figure 2(a). For this
case, where the only generalized stress different from zero is the bending moment,
the reference solution can be easily computed. In particular the free end vertical
displacement is:
W R2
f= , (14)
EI
where symbols which are not shown in Figure 2 are the Young’s modulus, E, and
the area moment of inertia of the cross-section, I.
Table 1 reports the free end vertical displacement by keeping fixed the number
of elements, ne = 8, and varying both the NURBS order p and the Gauss points
number ng in the range p, . . . , 4p. Reference value is obtained for a rectangular
cross-section with thickness b=0.2 m, depth h=0.01 m, and unitary values for W ,
R, E, i.e.: W =1 Nm, R=1 m, E=1 GPa.
These results suggest that a number of Gauss points ng equal to the NURBS
order p is the best choice, since it furnishes accurate results with a low computa-
tional cost. Finally, we noticed that when ng < p, i.e. the number of Gauss points
is less than the NURBS order, some ill-conditioning of the stiffness matrix occurs,
in particular for higher values of p. Therefore, for the sake of conciseness, only
values of ng ≥ p have been shown in Table 1. Moreover, unless otherwise noticed,
the choice ng = p has been adopted for all tests, whose results are shown in the
sequel.
Table 2 is again referred to the test reported in Figure 2(a) and reports the
numerical solution for the vertical displacement f of the free end; this time it
is either the number of elements, ne , or the NURBS order, p, which is going to
be increased. Besides the accuracy of the results, it is worth noticing that the
computational cost is essentially linked to the number of degrees of freedom (dofs)
of the problem. Since ne elements of order p give 3(ne + p) dofs, it is interesting
to compare results having the same computational cost, i.e. by keeping fixed the
6 ANTONIO CAZZANI, MARCELLO MALAGÙ, AND EMILIO TURCO

Table 1. Cantilever circular arch subjected to a couple at the free

end: vertical displacement of the free end by varying the NURBS
order p and the number of Gauss points ng (ne = 8, reference value
of the vertical displacement f = 6 cm).

ng
p p p+1 2p 4p
2 5.9975 4.7793 4.7806 4.7806
3 5.9996 5.9994 5.9994 5.9994
4 5.9999 5.9999 5.9999 5.9999
5 5.9999 5.9999 5.9999 5.9999

value of ne + p. We point out, to this end, that increasing the order p produces
better results than increasing ne .

Table 2. Cantilever circular arch subjected to a couple at the free

end: vertical displacement f of the free end for different number of
elements ne and NURBS order p (reference value f = 6 cm). The
adopted number of Gauss points is, in all shown cases, ng = p.

p
ne 2 3 4 5
4 5.9565 5.9852 5.9991 5.9999
8 5.9975 5.9996 5.9999 5.9999
16 5.9998 5.9999 5.9999 5.9999
32 5.9999 5.9999 5.9999 5.9999

Keeping in mind the observations reported in [15], we consider for the same test
the vertical displacement f as a function of the R/h ratio and the NURBS order
p, but keeping fixed the number of elements: in this case ne = 8, see Table 3. We
note that the results for low NURBS order p are poor when the R/h ratio increases
and this is a symptom of the occurrence of locking. To better investigate this
phenomenon, the behaviour of 2nd- and 3rd-order NURBS for the case R/h = 104
and R/h = 105 has been studied by increasing ne . From Table 4 we deduce that
also in this extreme cases, convergence is recovered by simply increasing ne .

Table 3. Cantilever circular arch subjected to a couple at the free

end: vertical displacement f of the free end, scaled by the reference
solution, varying R/h ratio and NURBS order p (ne = 8).

p
R/h 2 3 4 5
102 0.99958 0.99993 0.99998 0.99998
103 0.98815 0.99660 0.99995 0.99998
104 0.49402 0.86520 0,99702 0.99997
105 0.00968 0.06190 0.98163 0.99912
 ISOGEOMETRIC ANALYSIS OF PLANE–CURVED BEAMS 7

Table 4. Cantilever circular arch subjected to a couple at the free

end: vertical displacement f of the free, scaled by the reference
solution, varying ne for R/h = 104 . and R/h = 105 (p = 2 and
p = 3; ng = p).

p
R/h ne 2 3
16 0.98163 0.99500
104 32 0.99910 0.99985
64 0.99997 0.99998
16 0.37252 0.79427
105 32 0.97168 0.99198
64 0.99902 0.99838

As a second test we choose the ring under internal pressure depicted in Fig-
ure 2(b). In this case, the only non-zero generalized stress component is the axial
force N . Thanks to the polar symmetry it is straightforward to evaluate the refer-
ence solution. In particular, the axial force Nref and the radial displacement wref
are given by:
Nref = qR,
qR2 (15)
wref = .
EA
In Table 5 the radial displacement w and the axial force N are reported versus p
and ne (the number of elements being referred to a quarter of the whole ring). Both
values of w and N are scaled by the reference solution. The results are relative to
the following set of data: q = 1 kN/m, R = 1 m, E = 1 GPa and A = 1 cm2 .
First, also for this test the influence of the number of Gauss points has been
investigated. The results confirm that ng = p is again the best choice, even though
they are not reported here for space-saving reasons. Successively, the accuracy
of the results has been checked by computing N and w (properly scaled by the
reference solutions) as a function of the number of elements, ne , and of the NURBS
order, p.
Again, as in the previous test, the conclusion is that the effect of the so-called
p-refinement (i.e. increasing p) produces better results than the h-refinement (cor-
responding to increasing the number of elements, ne ). Indeed a relatively poor
mesh can provide accurate results as long as p is suitably increased.
We point out that convergence is not uniform: this comes out by the k-refinement
scheme implemented in our code for obtaining the numerical results. This kind of
refinement, which is typical of isogeometric analysis, is different from the h- and
p-refinement schemes of classical finite elements; see [1] and [18] for details.
Remembering the foregoing test and the loss of accuracy for small values of ne
and p, we also control the accuracy when the R/h ratio increases. Table 6 reports
the radial displacement w and the axial force N , both scaled for the respective
reference values, for the 2nd -order NURBS and ne = 4, when the R/h ratio varies
in the range 102 , . . . , 105 . In this case there is no loss of accuracy, even for high
values of R/h.
8 ANTONIO CAZZANI, MARCELLO MALAGÙ, AND EMILIO TURCO

Table 5. Ring under internal pressure: radial displacement w and

axial force N (both scaled for the reference solution) varying the
number of element ne and NURBS order p.

p
2 3 4 5
ne w/wref N/Nref w/wref N/Nref w/wref N/Nref w/wref N/Nref
2 0.99816 1.0035 1.0000 1.0000 0.99824 1.0000 1.0000 1.0000
4 0.99464 1.0012 0.99981 1.0000 0.99999 1.0000 1.0000 1.0000
8 0.99601 1.0002 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
16 0.99971 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
32 0.99998 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

Table 6. Ring under internal pressure: radial displacement w and

axial force N (both scaled for the reference solution) at varying
R/h values for NURBS order p = 2 (ne = 4).

R/h
102 103 104 105
w/wref 0.99464 0.99429 0.99429 0.99426
N/Nref 1.0012 1.0013 1.0013 1.0013

3.2. Cantilever circular arch under shear load at the free end. As a third
test we consider the cantilever circular arch, see Figure 3, loaded with a vertical unit
force P (in N) acting on the free end. We further assume the values R = 2 m for the
arc radius, E = 80 GPa for the Young’s modulus, ν = 0.2 for the Poisson’s ratio.
Finally, a rectangular cross-section with tickness b = 0.2 m and depth h = 0.01 m
is considered.

ψ
R

Figure 3. Sketch of the cantilever arch under shear load.

Also in this case, reference solution can easily be calculated and written as:
u = −P (c1 (sin ψ + ψ cos ψ) + c2 sin ψ), w = P c1 ψ sin ψ, ϕ = −P c3 sin ψ, (16)
where u, w, and ϕ are the tangential displacement, the radial displacement and the
cross-section rotation, respectively, expressed in terms of the angle ψ = s/R (see
 ISOGEOMETRIC ANALYSIS OF PLANE–CURVED BEAMS 9

again Figure 3) and of the following compliance coefficients:

 
1 R R R3 R R3 R2
c1 = + + , c2 = + , c3 = . (17)
2 EA GAT EI GAT EI EI
1
Generalized stresses N , T and M and strain energy Φ = 2 P w|π/2 , assume the
values:
1 2 π
N = −P cos ψ, T = P sin ψ, M = −P R cos ψ, Φ= P c1 . (18)
2 2
In Table 7 the vertical displacement f = w|π/2 of the free end as a function of the
NURBS order p and the number of elements ne is reported. Even in this case, if we
compare analyses performed with the same computational cost, i.e. with the same
value of ne + p, it comes out that the higher the degree p, the better the result is.

Table 7. Cantilever arch under shear load: vertical displacement

f of the free end as a function of the element number ne and
NURBS order p (reference value=4.7124 mm).

p
ne 2 3 4 5
1 4.4665 4.7059 4.7120 4.7124
2 3.8229 3,4257 4.6090 4.7080
4 3.8794 4.5343 4.7059 4.7122
8 4.6329 4.7076 4.7124 4.7124
16 4.7108 4.7124 4.7124 4.7124
32 4.7124 4.7124 4.7124 4.7124

Table 8 reports the only non-zero reactions at the clamped end: the vertical
force VA and bending moment MA for various NURBS order p and element num-
ber ne . Naturally, since our NURBS elements are based on a displacement-type
compatible formulation, it is expected that stress quantities exhibit relatively poor
performances; nonetheless good results can be obtained even with few elements by
suitably increasing the NURBS order p.
In order to compare the performances of NURBS to those of classical Lagrangian
interpolation, we report in Table 9 the vertical displacement of the free end f and
the reactions at the clamped end VA and MA for different order p and number
of element ne . We notice the superior computational performances of NURBS
interpolation with reference to the Lagrangian one.

3.3. Incomplete ring under vertical load. This numerical test was chosen since
it allows a comparison with a large number of finite elements reported in technical
literature. According to [19] and later [20], the problem sketched in Figure 4 was
analyzed using these data: load P = 1 lb; radius R = 2.935 in; Young’s modulus
E = 1.05 × 107 psi; Poisson’s ratio ν = 0.3; rectangular cross-section with thickness
b = 1.2 in and depth h = 0.125 in.
For this test we have also computed the analytical solution by integrating the
system of six linear first-order ODEs resulting from (1)–(3) and properly enforcing
the boundary conditions (BCs), account taken of symmetry with reference to the
10 ANTONIO CAZZANI, MARCELLO MALAGÙ, AND EMILIO TURCO

Table 8. Cantilever arch under shear load: reactions at the

clamped end vs. the element number ne and NURBS order p.

p ne VA [N] MA [Nm]
2 16 4.1624 2.0008
32 1.7861 2.0001
64 1.1960 2.0000
3 8 0.78097 1.9992
16 0.97405 1.9999
32 0.99680 2.0000
64 0.99960 2.0000
4 4 1.1809 1.9945
8 1.0044 2.0000
16 1.0001 2.0000
32 1.0000 2.0000
64 1.0000 2.0000
5 2 1.5499 1.9861
4 1.0041 2.0015
8 0.99903 2.0000
16 0.99997 2.0000
32 1.0000 2.0000
64 1.0000 2.0000
analytical 1 2

R
ψ
π/4

Figure 4. Sketch of incomplete ring under vertical load.

vertical axis. The generalized stresses N , T , M are given by (see again Figure 4
for the definition of the angle ψ):

N =A1 cos ψ + A2 sin ψ,

T =A2 cos ψ − A1 sin ψ, (19)
M =A2 R sin ψ − A1 R(1 − cos ψ) + A3 ,
 ISOGEOMETRIC ANALYSIS OF PLANE–CURVED BEAMS 11

Table 9. Cantilever arch under shear load: vertical displacement

f of the free end and reaction components VA and MA at the
clamped end as a function of the element number ne and of the
order of Lagrangian interpolation pL .

pL ne f [mm] VA [N] MA [Nm]

2 16 4.1337 320.22 1.6637
32 4.6685 96.007 1.9737
64 4.7098 25.068 1.9984
3 8 4.7087 0.83102 2.0124
16 4.7116 0.97569 2.0026
32 4.7122 0.99862 2.0005
64 4.7124 0.99994 2.0001
4 4 4.7124 8.5569 2.0116
8 4.7124 0.18327 2.0003
16 4.7124 0.94724 2.0000
32 4.7124 0.99669 2.0000
64 4.7124 0.99979 2.0000
5 2 4.7126 0.29883 1.9959
4 4.7124 1.0093 1.9999
8 4.7124 1.0003 2.0000
16 4.7124 1.0000 2.0000
32 4,7124 1.0000 2.0000
64 4.7124 1.0000 2.0000
analytical 4.7124 1 2

while the tangential and radial displacement and the cross-section rotation are:
u =A1 [c1 (ψ cos ψ + sin ψ) − c2 sin ψ − c3 R(sin ψ − ψ)]+
A2 [c1 ψ sin ψ − c3 R(1 − cos ψ)] + A3 c3 (sin ψ − ψ),
w =A2 [c1 (ψ cos ψ − sin ψ) + c2 sin ψ − c3 R sin ψ]− (20)
A1 [c1 ψ sin ψ − c3 R(1 − cos ψ)] − A3 c3 (1 − cos ψ),
ϕ =A1 c3 (sin ψ − ψ) + A2 c3 (1 − cos ψ) + A3 c3 ψ/R.
In (20) the same compliance coefficients c1 , c2 , c3 defined in (17) are used as a
short-hand notation, while constants A1 , A2 and A3 appearing in (19)–(20) assume
these values:
2P [c1 (α2 /R) − c3 (1 − cos α)] sin α
A1 = − ,
4c1 (α2 /R) − 2c3 (1 − cos 2α) + 2α sin 2α(c1 /R − c2 /R + c3 )
2P [c1 (α/R)(α cos α + sin α) − c2 (α/R) sin α − c3 sin α(sin α − α)]
A2 = , (21)
4c1 (α2 /R) − 2c3 (1 − cos 2α) + 2α sin 2α(c1 /R − c2 /R + c3 )
2P [c1 α(1 − cos α − α sin α) + (c1 − c2 ) sin α(cos α − 1)]
A3 = .
4c1 (α2 /R) − 2c3 (1 − cos 2α) + 2α sin 2α(c1 /R − c2 /R + c3 )
Setting α = 7π/8, for the given data the vertical displacement f of the incomplete
ring middle point results in f = 1.063161841 × 103 when all strain contributions are
taken into account; on the other hand, if shear deformation is discarded f |γ=0 =
1.059622233 × 103, while if both extensional and shear deformation are disregarded
12 ANTONIO CAZZANI, MARCELLO MALAGÙ, AND EMILIO TURCO

it results f |γ=0;=0 = 1.058121422 × 103 . The analytic solution provided by [19],

f = 1.06364 × 103 , does not seem compatible with the assumed geometric and
mechanical data.
Table 10 reports the vertical displacement f under the unitary load varying the
NURBS order p and the numbers of element ne on one-half of the structure. In the
same Table are reported both the analytical solution and the results obtained in [19]
and [20] with several type of finite elements. In order to simplify the comparison,
a column with the number of nodes nn is also reported. The Table shows the good
performances of the NURBS element family.

Table 10. Incomplete ring under vertical load: middle point ver-
tical displacement f .

p ne nn f × 10−3 Formulation Geometry

8 10 0.99727 compatibility exact
2 16 18 1.06201 ” ”
32 34 1.06314 ” ”
8 11 1.05863 ” ”
3 16 19 1.06311 ” ”
32 35 1.06316 ” ”
8 12 1.06301 ” ”
4 16 20 1.06316 ” ”
32 36 1.06316 ” ”
8 13 1.06315 ” ”
5 16 21 1.06316 ” ”
32 37 1.06316 ” ”
analytical — — 1.06316 — —
[19] 8 16 1.05641 compatibility exact
[19] 8 16 1.05641 hybrid exact
[19] 8 16 1.05619 mixed exact
[20] 1 3 1.07069 equilibrium cubic parametric
[20] 2 6 1.06384 equilibrium cubic parametric

3.4. Clamped-clamped semi-circular arch under distributed load. Now,

we consider the clamped-clamped semi-circular arch depicted in Figure 5. This
time, we assume R = 1 m, a square cross-section with side b = 0.1 m, Young’s
modulus E = 1 GPa, Poisson’s ratio ν = 0 and, finally, a load q = 1 N/m. The
analytical solution is provided again by integrating the system (1)–(3), properly
enforcing BCs and symmetry with reference to the vertical axis. In this case, by
taking into account that the load is uniformly distributed per unit of projection,
it results: qt = −q sin ψ cos ψ and qr = q sin2 ψ. The generalized stresses N , T , M
and displacements u, w and ϕ are, respectively, given by:

N =A2 sin ψ − qR cos2 ψ,

T =A2 cos ψ + qR cos ψ sin ψ, (22)
2
M =A2 R sin ψ + A3 − qR /2(1 + 1/2 cos 2ψ),
 ISOGEOMETRIC ANALYSIS OF PLANE–CURVED BEAMS 13

and
u =A2 [c1 ψ sin ψ − c3 R(1 − cos ψ)] − A3 c3 (ψ − sin ψ)+
A5 sin ψ − qR[sin 2ψ(2/3c1 − 1/6c2 − 1/8c3 R) − ψc3 R/2],
w =A2 [c1 (ψ cos ψ − sin ψ) + c2 sin ψ − c3 R sin ψ] − A3 c3 (1 − cos ψ)+ (23)
A5 cos ψ + qR[c1 − 1/2c2 + 1/2c3 R − cos 2ψ(1/3c1 + 1/6c2 − 1/4c3 R)],
ϕ =A3 ψ/Rc3 − A2 c3 (cos ψ − 1) − qRc3 (ψ/2 + 1/8 sin 2ψ).
As before, in (23) the compliance coefficients c1 , c2 , c3 defined in (17) have been
used for brevity, while constants A2 , A3 and A5 appearing in (22)–(23) assume the
values:
8πq(c1 − c2 ) + 3πqRc3
A2 = ,
6π 2 (c1 /R) − 24c3
qR2 16πqR(c1 − c2 ) + 6πqR2 c3
A3 = − , (24)
2 6π 3 (c1 /R) − 24πc3
2qR(c1 − c2 ) 3qR2 c3
A5 =− − .
3 4
For the given data, the vertical displacement of the arch middle point, evaluated
in exact form with a Computer Algebra System is vC = 1.018188371 mm when
all strain contribution are taken into account; if shear deformation is discarded
vC |γ=0 = 0.982381687 mm, while if both extensional and shear deformation are
disregarded it results vC |γ=0;=0 = 0.817230149 mm. The reactions at point A are
instead: HA = 0.554438 N, VA = 1.000000 N and MA = 0.102966 Nm.
We point out that some results of the analytical solution provided by [20], i.e.
vC = 1.02086 mm, HA = 0.554915 N, MA = 0.103980 Nm differ from the analytical
solution presented above.

R
A ψ

Figure 5. Sketch of the clamped-clamped semi-circular arch un-

der distributed load.

Table 11 reports the values of the vertical displacement vC of the point C and
the reactions at the clamped end A: the horizontal, HA , and vertical, VA , force
components and bending moment, MA , for different NURBS order p and number
of elements ne (referred to one half of the arch), along with our analytical values. In
the same table, the results obtained by [20] neglecting transverse shear deformation
are reported.
We remark the high accuracy of our results, for both displacements and gener-
alized stresses even if convergence rate for stresses is slower and a large number of
14 ANTONIO CAZZANI, MARCELLO MALAGÙ, AND EMILIO TURCO

elements or a higher NURBS degree is required to achieve a good quality approxi-

mation.

Table 11. Clamped-clamped semi-circular arch under distributed

load: middle vertical displacement and clamp reactions.

p ne vC [mm] HA [N] VA [N] MA [Nm]

2 16 1.01812 0.728620 0.932677 0.102130
32 1.01819 0.598888 0.981880 0.102778
3 8 1.01811 0.571273 0.968679 0.103210
16 1.01819 0.556235 0.996308 0.102981
32 1.01819 0.554667 0.995452 0.102968
4 8 1.01819 0.553398 1.00105 0.102945
16 1.01819 0.554405 1.00003 0.102966
32 1.01819 0.554437 1.00000 0.102966
5 8 1.01819 0.554595 0.999814 0.102966
16 1.01819 0.554440 0.999995 0.102966
32 1.01819 0.554438 1.000000 0.102966
analytical — 1.01819 0.554438 1 0.102966
[20] 1 1.02135 0.555009 1.0 0.104159
[20] 2 9.9091 0.500476 1.0 0.092636
[20] 4 1.00177 0.540775 1.0 0.101238

3.5. Three-hinged lancet arch under self weight. The three-hinged pointed
arch reported in Figure 6 is taken into consideration with the following set of data:
E = 1 GPa, ν = 0.2, R = 1 m, b = 0.2 m, h = 0.01 m and q0 = 1 kN/m.

C
q0
R

A B π/4
R ψ

Figure 6. Sketch of the three-hinged lancet arch under self weight.

As before, the analytical solution has been computed by integrating (1)–(3) and
imposing the relative boundary and symmetry conditions. In particular at the arch
tip, i.e. for α = π/4, the following BCs, which involve the bending moment M ,
 ISOGEOMETRIC ANALYSIS OF PLANE–CURVED BEAMS 15

the horizontal displacement h, and the vertical component of the internal force V ,
must be met:
M |α = 0; h|α = u|α sin α + w|α cos α = 0; V |α = N |α cos α − T |α sin α = 0.
Since the load is uniformly distributed per unit of arc-length (and not per unit
of projection), it results qt = −q0 cos ψ and qr = q0 sin ψ. Then, the generalized
stresses N , T , M are given by:
N = A1 cos ψ + A2 sin ψ + q0 Rψ cos ψ,
T = A2 cos ψ − A1 sin ψ − q0 Rψ sin ψ, (25)
2
M = A1 R(cos ψ − 1) + A2 R sin ψ + q0 R (ψ cos ψ − sin ψ),
while displacements and cross-section rotation are:
u = A1 [c1 ψ cos ψ + (c1 − c2 − c3 R) sin ψ + c3 Rψ] + A4 cos ψ+
A2 [c1 ψ sin ψ + c3 R(cos ψ − 1)] + A6 R(cos ψ − 1)+
q0 [c1 R(ψ 2 /2) cos ψ + (c1 − c2 )R(cos ψ + 2ψ sin ψ)/4−
c3 R2 (2 cos ψ + ψ sin ψ)],
w = A1 [−c1 ψ sin ψ + c3 R(1 − cos ψ)] − A4 sin ψ+ (26)
A2 [c1 (ψ cos ψ − sin ψ) + c2 sin ψ − c3 R sin ψ] − A6 R sin ψ+
q0 [(c2 − c1 )R(2ψ cos ψ − sin ψ)/4 − c1 R(ψ 2 /2) sin ψ+
c3 R2 (sin ψ − ψ cos ψ)],
ϕ = A1 c3 (sin ψ − ψ) + A2 c3 (1 − cos ψ) + A6 + q0 c3 R(2 cos ψ + ψ sin ψ).
In (26) the compliance coefficients c1 , c2 , c3 defined in (17) have been used to
simplify notation, while constants A1 , A2 , A4 and A6 appearing in (25)–(26) assume
these values:
A1 = − q0 Rα,
A2 = q0 R(1 − α/ sin α),
A4 = q0 R[2c3 R + (c2 − c1 )/4], (27)
2 2
A6 = − q0 [c1 α / sin α + (c1 − c2 )(1/2 cos α − α/ tan α)−
(c1 + c2 )α/(2 sin α) + c3 R(1 + α2 + α/ tan α − α/ sin α + cos α)].
For the given data, the vertical displacement of the arch middle point, evaluated
in exact form with a Computer Algebra System is vC = 5.47802398 mm when all
strain contribution are taken into account. The reactions at point A are instead:
HA = 110.72073454 N, VA = 785.39816340 N.
Table 12 reports the values of the vertical displacement vC and the reactions at
support A, the horizontal, HA , and the vertical, VA , force components for different
NURBS order p and number of element ne (only one half of the arch was considered
due to the symmetry).

4. Concluding remarks
In this paper, NURBS finite elements for the analysis of plane curved beams
have been investigated. Numerical results show, for this family of finite elements,
good performances. In particular, accuracy and convergence ratio by varying both
the number of elements and the NURBS order have been thoroughly studied. The
16 ANTONIO CAZZANI, MARCELLO MALAGÙ, AND EMILIO TURCO

Table 12. Three-hinged lancet arch under self weight: tip dis-
placement and reaction forces at the support A.

p ne vC [mm] HA [kN] VA [kN]

2 32 5.47801 0.314991 0.790041
3 16 5.47800 0.146215 0.771868
32 5.47802 0.114962 0.783706
8 5.47802 0.115145 0.782591
4 16 5.47802 0.110882 0.785265
32 5.47802 0.110727 0.785391
8 5.47802 0.110620 0.785651
5 16 5.47802 0.110723 0.785402
32 5.47802 0.110721 0.785398
analytical — 5.47802 0.110721 0.785398

influence of the number of Gauss points on the results was considered too, and
a simple rule was proposed, which consists of using p Gauss points to integrate
the stiffness matrix corresponding to a p-th order NURBS finite element. Besides
numerical results, some new analytical solutions have been provided in order to
expand those available in the technical literature. This paper has surely given
some useful guidelines about the isogeometric approach for plane curved beams
but there are a few aspects which, in our opinion, deserve further extensions. What
we consider the most interesting ones are listed below:
(1) The first extension descends from observing the generalized stress results
produced by the numerical model. Since they are clearly less accurate than
the corresponding displacements, a possible improvement could be found
in the framework of mixed formulations, like those proposed in [21], or in
special stress recovery techniques (see [22]).
(2) NURBS interpolation is somewhat different from polynomial interpolation,
so that using the same integration rules, (i.e. standard Gauss formulae
which rigorously apply only to polynomials) might not be the best way
to obtain accurate results or the most suitable from the point of view of
computational cost. Different solutions have been proposed, for example,
in [23]. For these reasons, we think that this issue would require additional
investigations.
(3) Some possible applications of the Timoshenko beam model and non-classical
applications of 1D beam-like model are worth mentioning. Among them
there is the analysis of beams made of composite and functionally graded
materials where proper relations for 1D elastic properties are established
([24]); the extension of the Timoshenko beam model to the thermal prob-
lems for beams within the framework of two temperature model ([25]); the
description of particular phenomena like a line tension of phase interfaces
in 2D structures such as a martensitic film ([26]).
(4) The promising results obtained for plane curved beams suggest an extension
to space-curved beams and also to more complicated two-dimensional struc-
tures such as shells. Moreover, extensions to investigate problems where
both geometric and material nonlinearities can occur ([27]), or to linear and
non-linear structural dynamics ([28], [29]) and buckling ([30], [31]) could be
 ISOGEOMETRIC ANALYSIS OF PLANE–CURVED BEAMS 17

achieved, as well as the search for fundamental solutions through Green’s

operator ([32]).
(5) Some applications require considering continuum models which are richer
than the standard Cauchy one: it would be possible then trying to extend
the proposed isogeometric approach to second and higher-order continua
(see [33], [34], [35], [36]).
(6) A foreseable application is the study of larger systems made of many
NURBS special elements, for example for the analysis of damage mechanics
in the framework of structural health monitoring (see [37], [38], [39], [40])
or along the statistical approach developed for Hookean spring networks
(see [41], [42], [43]).
(7) Finally, an interesting application would involve tackling composite or porous
([44]) or fibrous materials ([45]), even in cases when complex physical
problems like bio-medical ones whose evolution is governed by integro-
differential operators and material properties are subject to sudden changes
([46], [47], [48]). In all these cases it is believed that the ability of properly
modulating the regularity of the basis function through the use of NURBS
could make more tractable such problems.

Acknowledgements. The financial support of RAS, the Autonomous Region of

Sardinia, under grant number F71J09999350002-CRP1 475 (Legge Regionale 7/2007,
bando 2008, Progetto MISC: Metodi Isogeometrici per Strutture Curve) is gratefully
acknowledged.

References
[1] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs. Isogeometric analysis: CAD, finite elements,
NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and
Engineering, 194:4135–4195, 2005.
[2] J. A. Cottrell, T. J. R. Hughes, and Y. Bazilevs. Isogeometric Analysis: Toward Integration
of CAD and FEA. Wiley, 2009.
[3] T. J. R. Hughes, A. Reali, and G. Sangalli. Efficient quadrature for NURBS-based isoge-
ometric analysis. Computer Methods in Applied Mechanics and Engineering, 199:301–313,
2010.
[4] D. J. Benson, Y. Bazilevs, M. C. Hsu, and T. J. R. Hughes. Isogeometric shell analysis: the
Reissner-Mindlin shell. Computer Methods in Applied Mechanics and Engineering, 199:276–
289, 2010.
[5] A. Benedetti and A. M. Tralli. A new hybrid F.E. model for arbitrarily curved beam–I. linear
analysis. Computer and Structures, 33(6):1437–1449, 1989.
[6] C. Gontier and C. Vollmer. A large displacement analysis of a beam using a CAD geometric
definition. Computer and Structures, 57(6):981–989, 1995.
[7] R. Echter and M. Bischoff. Numerical efficiency, locking and unlocking of NURBS finite
elements. Computer Methods in Applied Mechanics and Engineering, 199(5–8):374–382, 2010.
[8] L. Beirão da Veiga, C. Lovadina, and A. Reali. Avoiding shear locking for the Timoshenko
beam problem via isogeometric collocation methods. Computer Methods in Applied Mechan-
ics and Engineering, 241–244:38–51, 2012.
[9] F. Auricchio, L. Beirão da Veiga, J. Kiendl, C. Lovadina, and A. Reali. Locking-free iso-
geometric collocation methods for spatial Timoshenko rods. Computer Methods in Applied
Mechanics and Engineering, 263:113–126, 2013.
[10] M. Cuomo and L. Greco. Isogeometric analysis of space rods: considerations on stress locking.
In ECCOMAS 2012, pages 1–19, 2012.
[11] L. Greco and M. Cuomo. B-Spline interpolation of Kirchhoff-Love space rods. Computer
Methods in Applied Mechanics and Engineering, 256:251–269, 2013.
18 ANTONIO CAZZANI, MARCELLO MALAGÙ, AND EMILIO TURCO

[12] L. Greco and M. Cuomo. A locking-free multipatch B-spline element for the analysis of curved
3D rod elements. In Aimeta 2013, pages 1–10, 2013.
[13] M. Cuomo, L. Contrafatto, and L. Greco. A variational model based on isogeometric inter-
polation for the analysis of cracked bodies. International Journal of Engineering Sciences,
page (in press), 2013.
[14] L. Greco and M. Cuomo. An implicit G1 multi patch B-spline interpolation for Kirchhoff-
Love space rods . Computer Methods in Applied Mechanics and Engineering, 269:173–197,
2014.
[15] R. Bouclier, T. Elguedj, and A. Combescure. Locking free isogeometric formulations of curved
thick beams. Computer Methods in Applied Mechanics and Engineering, 245–246:144–162,
2012.
[16] L. Piegl and W. Tiller. The NURBS Book. Springer, 2nd edition, 1996.
[17] C. de Falco, A. Reali, and R. Vazquez. GeoPDEs: a research tool for isogeometric analysis
of PDEs. Advances in Engineering Software, 42(12):1020–1034, 2011.
[18] J. A. Cottrell, T. J. R. Hughes, and A. Reali. Studied of refinement and continuity in iso-
geometric structural analysis. Computer Methods in Applied Mechanics and Engineering,
196:4160–4183, 2007.
[19] H. Stolarski and T. Belytschko. Shear and membrane locking in curved C0 elements. Com-
puter Methods in Applied Mechanics and Engineering, 41(3):279–296, 1983.
[20] L. Molari and F. Ubertini. A flexibility-based element for linear analysis of arbitrarily curved
arches. International Journal for Numerical Methods in Engineering, 65:1333–1353, 2006.
[21] M. Cuomo and G. Ventura. Complementary energy approach to contact problems based
on consistent augmented Lagrangian formulation. Mathematical and Computer Modelling,
28:185–204, 1998.
[22] D. Ciancio, I. Carol, and M. Cuomo. On inter-element forces in the FEM-displacement for-
mulation, and implications for stress recovery. International Journal for Numerical Methods
in Engineering, 66:502–528, 2006.
[23] M. J. Borden, M. A. Scott, J. A. Evans, and T. J. R. Hughes. Isogeometric finite element
data structures based on Bézier extraction of NURBS. International Journal for Numerical
Methods in Engineering, 87(1–5):15–47, 2011.
[24] M. Bı̂rsan, H. Altenbach, T. Sadowski, V.A. Eremeyev, and D. Pietras. Deformation anal-
ysis of functionally graded beams by the direct approach. Composites Part B: Engineering,
43:1315–1328, 2012.
[25] Holm Altenbach, Mircea Bı̂rsan, and Victor A. Eremeyev. On a thermodynamic theory of
rods with two temperature fields. Acta Mechanica, 223:1583–1596, 2012.
[26] W. Pietraszkiewicz, V. Eremeyev, and V. Konopińska. Extended non-linear relations of elastic
shells undergoing phase transitions. ZAMM—Zeitschrift für Angewandte Mathematik und
Mechanik, 87:150–159, 2007.
[27] G. Oliveto and M. Cuomo. Incremental analysis of plane frames with geometric and material
nonlinearities. Engineering Structures, 10:2–12, 1988.
[28] A. Luongo, G. Rega, and F. Vestroni. On nonlinear dynamics of planar shear indeformable
beams. Journal of Applied Mechanics, Transactions ASME, 53:619–624, 1986.
[29] A. Di Egidio, A. Luongo, and A. Paolone. Linear and non-linear interactions between static
and dynamic bifurcations of damped planar beams. International Journal of Non-Linear
Mechanics, 42:88–98, 2007.
[30] A. Luongo. On the amplitude modulation and localization phenomena in interactive buckling
problems. International Journal of Solids and Structures, 27:1943–1954, 1991.
[31] A. Luongo. Mode localization in dynamics and buckling of linear imperfect continuous struc-
tures. Nonlinear Dynamics, 25:133–156, 2001.
[32] M. Cuomo and G. Ventura. An explicit formulation of the Green’s operator for general one-
dimensional structures. European Journal of Mechanics, A/Solids, 21:493–512, 2002.
[33] F. dell’Isola, G.C. Ruta, and R.C. Batra. Second-order solution of Saint-Venant’s problem
for an elastic pretwisted bar using Signorini’s perturbation method. Journal of Elasticity,
49:113–127, 1997.
[34] R.C. Batra, F. dell’Isola, and S. Vidoli. A second-order solution of Saint-Venant’s problem
for a piezoelectric circular bar using Signorini’s perturbation method. Journal of Elasticity,
52:75–90, 1998.
 ISOGEOMETRIC ANALYSIS OF PLANE–CURVED BEAMS 19

[35] J.-J. Alibert, P. Seppecher, and F. dell’Isola. Truss modular beams with deformation energy
depending on higher displacement gradients. Mathematics and Mechanics of Solids, 8:51–73,
2003.
[36] F. dell’Isola, P. Seppecher, and A. Madeo. How contact interactions may depend on the shape
of Cauchy cuts in N-th gradient continua: Approach à la D’Alembert. ZAMP—Zeitschrift
für Angewandte Mathematik und Physik, 63:1119–1141, 2012.
[37] N. Roveri and A. Carcaterra. Damage detection in structures under traveling loads by Hilbert-
Huang transform. Mechanical Systems and Signal Processing, 28:128–144, 2012.
[38] N. Roveri and A. Carcaterra. Structural health monitoring of time-varying systems by output-
only identification. In Proceedings of ISMA 2012, International Conference on Noise and
Vibration Engineering, Leuven, September 2012. KU Leuven. ISBN: 9789073802896.
[39] N. Roveri, A. Carcaterra, and A. Sestieri. Integrated system for shm and wear estimation of
railway infrastructures. In Surveillance 7 — 7th International Conference on Acoustical and
Vibratory Surveillance Methods and Diagnostic Techniques, pages P39/1–19, Chartres, 2013.
[40] A. Sestieri, A. Carcaterra, and N. Roveri. New monitoring technologies in mechanical systems.
In Surveillance 7 — 7th International Conference on Acoustical and Vibratory Surveillance
Methods and Diagnostic Techniques, pages P4/1–21, Chartres, 2013.
[41] A. Rinaldi and Y.-C. Lai. Statistical damage theory of 2D lattices: Energetics and physical
foundations of damage parameter. International Journal of Plasticity, 23:1796–1825, 2007.
[42] A. Rinaldi. Rational damage model of 2D disordered brittle lattices under uniaxial loadings.
International Journal of Damage Mechanics, 18(3):233–257, 2009.
[43] A. Rinaldi and L. Placidi. A microscale second gradient approximation of the damage param-
eter of quasi-brittle heterogeneous lattices. ZAMM—Zeitschrift für Angewandte Mathematik
und Mechanik, pages 1–16, 2013. DOI:10.1002/zamm.201300028.
[44] A. Madeo, F. dell’Isola, and F. Darve. A continuum model for deformable, second gradient
porous media partially saturated with compressible fluids. Journal of the Mechanics and
Physics of Solids, 61:2196–2211, 11 2013.
[45] M. Ferretti, A. Madeo, F. dell’Isola, and P. Boisse. Modeling the onset of shear boundary
layers in fibrous composite reinforcements by second-gradient theory. ZAMP—Zeitschrift für
Angewandte Mathematik und Physik, pages 1–26, 2013. DOI:10.1007/s00033-013-0347-8.
[46] A. Madeo, T. Lekszycki, and F. dell’Isola. A continuum model for the bio-mechanical in-
teractions between living tissue and bio-resorbable graft after bone reconstructive surgery.
Comptes Rendus Mécanique, 339:625–640, 2011.
[47] T. Lekszycki and F. dell’Isola. A mixture model with evolving mass densities for describing
synthesis and resorption phenomena in bones reconstructed with bio-resorbable materials.
ZAMM—Zeitschrift für Angewandte Mathematik und Mechanik, 92:426–444, 2012.
[48] A. Madeo, D. George, T. Lekszycki, M. Nierenberger, and Y. Rémond. A second gradient
continuum model accounting for some effects of micro-structure on reconstructed bone re-
modelling. Comptes Rendus Mécanique, 340:575–589, 2012.

(Corresponding author) Dipartimento di Ingegneria Civile, Ambientale e Architettura

(DICAAR), Università degli Studi di Cagliari, Italia
E-mail address: [email protected]

Dipartimento di Ingegneria Civile, Ambientale e Architettura (DICAAR), Università

degli Studi di Cagliari, Italia
E-mail address: [email protected]

Dipartimento di Architettura, Design e Urbanistica (DADU), Università degli Studi

di Sassari, Italia
E-mail address: [email protected] / [email protected]