# Spatial Finite Element Analysis for Dynamic Response of Curved Thin-Walled Box Girder Bridges

Authors: Yinhui Wang; Yidong Xu; Zheng Luo; Haijun Wu; Liangliang Yan
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

According to the flexural and torsional characteristics of curved thin-walled box girder with the effect of initial curvature, 7 basic displacements of curved box girder are determined. And then the strain-displacement calculation correlations were established. Under the curvilinear coordinate system, a three-noded curved girder finite element which has 7 degrees of freedom per node for the vibration characteristic and dynamic response analysis of curved box girder is constructed. The shape functions are used as the interpolation functions of variable curvature and variable height to accommodate to the variation of curvature and section height. A MATLAB numerical analysis program has been implemented.


---

## Body

## 1. Introduction

The high torsional strength of box girder makes it particularly suitable for sharp curve alignment so as to accommodate to topography along traffic lines, skewed piers and abutment, superelevation, and transitions such as interchange ramp structures. Such accommodation plus aesthetic appearance makes it popular in curved bridges.Currently the static response of curved thin-walled box girder has been studied thoroughly. Comparatively speaking, the corresponding dynamic response analysis is more complicated due to the spatial coupling effect between bending and torsion and so there are still many potential as well as crucial problems regarding dynamic response, for example, forced vibration due to moving loads and earthquake, vehicle-bridge coupling vibration, and wind-induced vibration, for us to explore in depth [1–3].As for the numerical methods, dynamic finite element analysis of curved girder was developed on static finite element analysis basis. Much research has been conducted about the spatial mechanical behavior of curved box girder utilizing one-dimensional (denoted by 1D hereafter) curved box beam element. Y. Kim and Y. Y. Kim [4] carried out a higher-order 1D analysis of the in-plane flexural deformation of thin-walled curved box beams by adding an additional degree of freedom accounting for in-plane distortional deformation of cross section. Further, Y. Y. Kim and Y. Kim [5] proposed a five degrees of freedom (torsion, out-of-plane, flexure, warping, and distortion) 1D theory for the analysis of thin-walled curved rectangular box beams subjected to torsion and out-of-plane bending by including two additional degrees corresponding to warping and distortion. Cheng-Long et al. [6] put forward 1D finite element of 14 degrees of freedom at every node under convected cylinder coordinates for the shear lag effect analysis, which takes account of tension, compression, bending, torsion, warping, shear lag effect, and their coupling effect. Wang and Chen [7] constructed an iso-parameter finite element of three nodes and nine degrees of freedom per node to perform spatial analysis for box girder of variable curvature with initial curvature. Tension, compression, bending, torsion, warping, and distortion were all included. It applies to spatial static analysis as well as time-effect analysis and has many advantages over solid element and plate/shell element such as simpler model, less computation, more intuitive result representation, and wider application whether for dynamic property analysis or dynamic response analysis. Calik-Karakose et al. [8] constructed a four-noded curved shell finite element, which have 7 degrees of freedom per node, for the geometrically nonlinear analysis of beams curved in plan. The curved box girder structure was conceived as a sequence of macroelements having the form of transversal segments of identical topology where each slice is formed using a number of the curved shell elements. Using thin-walled elements with both warping torsion and distortion taken into account, Huang et al. [9] analyzed vibration of horizontally curved box girder bridges due to vehicle. Taking account of the effect of shear deformation and that of rotary inertias due to bending and torsional vibration, Wu and Chiang [10] derived the stiffness matrix and mass matrix of curved beam based on the local coordinate system and then applied them to analyze the natural vibration characteristics and dynamic response due to a moving load. N.-I. Kim and M.-Y. Kim [11] evaluated an exact dynamic element stiffness matrix of curved beam for the spatially coupled free vibration analysis of shear deformable thin-walled nonsymmetric curved beam subjected to initial axial force.Extending the spatial curved beam finite element theory considering effect of initial curvature developed earlier [7], the objective of this paper is to propose a 1D curved box girder dynamic finite element of three nodes and twenty-one degrees of freedom taking account of tension, compression, bending, torsion, and warping and then apply it to analyze the natural vibration characteristics and dynamic response. According to the flexural and torsional characteristics of curved thin-walled box girder with the effect of initial curvature taken into account and distortion ignored, seven basic displacements of curved box girder are determined. And then the relations between such basic displacements and displacements at generic point within the cross section were established. Furthermore, the strain (including normal strain and shear strain) displacement relationship formulas are deduced. Under the curvilinear coordinate system, a curved girder finite element with 3 nodes which has 7 degrees of freedom per node has been constructed and the corresponding shape functions have been constructed. The expression of element kinetic energy in terms of the first-order derivatives of the basic displacements with respect to time t was deduced, from which the consistent mass matrix was obtained. Afterwards, Rayleigh damping matrix, which is a liner combination of stiffness matrix and mass matrix, is introduced. Solving the frequency equations of undamped free vibration gave the natural frequencies and modes of structure. And then the dynamic responses of structure are worked out utilizing modal superposition method. In order to accommodate to the variation of curvature of element and the height of cross section, the shape functions are also used as the interpolation functions of variable curvature and variable-height. A MATLAB procedure based on the finite element formulation established herein is developed to implement numerical analysis. Finally, two case studies, one is for straight box girder and the other is for curved box girder, are conducted to confirm the validity and accuracy of the finite element put forward herein.


## 2. The Coordinate System and Basic Assumptions

The object to be studied in this paper is a curved box girder of variable curvature and variable height of cross section, whose cross section is trapezoidal and symmetric vertically. To address the analysis hereafter, we set up a curvilinear orthogonal coordinate system shown as Figure1, in which the positive directions of x-axis and y-axis are taken as the directions directing towards the center of the circle and downwards vertically, respectively. The horizontal curve along the centerline of the surface of top flange is taken as z-axis. The curvature (denoted by κ) and height (denoted by h) are variable along z-axis. Under such a coordinate, y
S and y
M are the ordinates of centroid of cross section (denoted by S) and torsional center (denoted by M), respectively. Due to initial curvature, the centroid of cross section does not coincide with the flexural center at e
,
y
S
,
z and herein e can be obtained by referring to the literature [12].Figure 1
Cross section and curvilinear orthogonal coordinate system of curved box girder.To facilitate our discussion, some assumptions based on which the analysis is performed should be introduced first as follows.(1)
The following factors should be taken into account.
(a)
The misalignment between cross section centroid andflexural center due to initial curvature of curved box girder.
(b)
The influence of shear strain due to deflection on elastic strain energy (the theory of Timoshenko beam).
(c)
The difference of deflection angle at every point of cross section due to initial curvature.
(2)
It is hypothesized that the normal stress and shear stress of cross section due to warping is uniform along the thickness of thin wall and the shear deformation of curved surface of open cross section is neglected.
(3)
The distribution of longitudinal displacement due to warping is assumed to be the same as that due to free torsion. An independent displacement variable due to warping,β
(
z
), is introduced.
(4)
Generally, the slabs of reinforced concrete box girder are thick and there are haunches to resist distortion. The diaphragm plates of steel box girder have the same function. Therefore, the influence of distortion of curved box girder on dynamic response is negligible enough to be neglected.


## 3. The Displacement and Strain at Generic Point within Curved Box Girder

### 3.1. The Displacement at Generic Point within Curved Box Girder

Denote the displacements at the flexural center parallel tox-axis, y-axis, and z-axis by u
c
(
z
), v
c
(
z
), and w
c
(
z
), respectively. Denote the deflection angles about radial neutral axis and vertical neutral axis of cross section by φ
x
(
z
) and φ
y
(
z
), respectively. Denote the torsion angle about torsional center and warping displacement by θ
(
z
) and β
(
z
).Considering the influence of shear strain due to deflection, we can expressφ
x
(
z
) and φ
y
(
z
) as(1)
φ
x
z
=
-
∂
v
c
∂
z
+
γ
x
φ
y
z
=
∂
u
c
∂
z
+
γ
y
,where γ
x and γ
y are the shear strains of y-direction and x-direction due to deflection, respectively.When the curvature is great,φ
x also varies with x and thus it should be written as(2)
φ
x
x
,
z
=
-
∂
v
c
∂
z
1
+
γ
x
.It follows from∂
z
1
=
(
R
0
-
x
/
R
0
)
∂
z that(3)
φ
x
x
,
z
≈
R
0
R
0
-
x
φ
x
z
=
1
+
κ
x
φ
x
z
,which reflects the influence of initial curvature onφ
x and where κ is the curvature of radial direction satisfying κ
=
1
/
R
0.And then, neglecting distortion, we can express the displacement at generic point(
x
,
y
,
z
) within curved box girder in terms of such seven basic displacements, say, u
c
(
z
), v
c
(
z
), w
c
(
z
), φ
x, φ
y, θ
(
z
), and β
(
z
) as (4)
u
x
,
y
,
z
=
u
c
z
-
y
-
y
M
θ
z
v
x
,
y
,
z
=
v
c
z
+
x
θ
z
w
x
,
y
,
z
=
w
c
z
-
φ
y
z
+
κ
w
c
z
x
-
e
+
y
-
y
S
1
+
κ
x
φ
x
z
-
ω
β
′
z
-
κ
φ
x
z
,where e is the eccentricity of flexural center with respect to centroid in the x-direction; ω is the principle sectorial coordinate of the cross section box girder; y
S and y
M are the coordinates of centroid and shear center in y-direction, respectively.


### 3.2. The Relation between the Longitudinal Normal Strain inz-Direction at Generic Point within the Box Girder and Seven Basic Displacements

LettingR
0 and R
0
-
x to be the longitudinal radius of curvature at the centroid (i.e., x
=
0) and generic point of the cross section of box girder, respectively, we can obtain(5)
∂
z
1
=
R
0
-
x
∂
z
R
0
,where ∂
z
1 and ∂
z are the partial derivative of z
1 and z.The abscissa isx, respectively. And herein z
1-axis is parallel to z-axis. Therefore, the longitudinal normal strain in z-direction at generic point within curved box girder considering the effect of initial curvature can be written as [7](6)
ε
z
=
R
0
R
0
-
x
·
∂
w
∂
z
-
u
R
0
-
x
=
1
+
κ
x
∂
w
∂
z
-
κ
u
.Substituting (4) into (6) yields(7)
ε
z
=
1
+
κ
x
w
C
′
-
φ
y
+
κ
w
c
x
-
e
′
+
1
+
κ
x
y
-
y
S
φ
x
′
-
ω
β
′
-
κ
φ
x
′
-
κ
u
c
-
y
-
y
M
θ
.The principle sectorial coordinateω in (7) should be discussed for such cases as top deck, bottom deck, right web, left web, right flange, and left flange, respectively. Due to the effect of variable curvature and variable height, h, κ, y
S, y
M, e, and ω are all functions of z.


### 3.3. The Relations between the Shear Strain at Generic Point within Curved Box Girder and Seven Basic Displacements

Considering any infinitesimal curved element ofd
z by d
s in top flange, bottom flange, or web, we can obtain the shear strain (
γ
) between the plane s
z with curvature taken into account due to a deformation w in the z-direction and u
s:(8)
γ
=
∂
w
∂
s
+
1
+
κ
x
∂
u
s
∂
z
,where u
s (see Figure 2) is the displacement measured along the center line of the closed cross section denoted by s. For the cases of top flange, right web, bottom flange, and left web, u
s is given as (9)
top  flange:
u
s
=
u
x
,
t
1
2
,
z
right  web:
u
s
=
-
u
x
,
y
,
z
cos
⁡
α
+
v
x
,
y
,
z
sin
⁡
α
bottom  flange:
u
s
=
-
u
x
,
h
+
t
1
2
,
z
left  web:
u
s
=
-
u
x
,
y
,
z
cos
⁡
α
-
v
x
,
y
,
z
sin
⁡
α
,where α is the inclination of web.Figure 2
u
s of curved cross section.Substituting (9) into (8) yields the relations of shear stain of top flange, right web, bottom flange, and left web with the displacements, respectively.


## 3.1. The Displacement at Generic Point within Curved Box Girder

Denote the displacements at the flexural center parallel tox-axis, y-axis, and z-axis by u
c
(
z
), v
c
(
z
), and w
c
(
z
), respectively. Denote the deflection angles about radial neutral axis and vertical neutral axis of cross section by φ
x
(
z
) and φ
y
(
z
), respectively. Denote the torsion angle about torsional center and warping displacement by θ
(
z
) and β
(
z
).Considering the influence of shear strain due to deflection, we can expressφ
x
(
z
) and φ
y
(
z
) as(1)
φ
x
z
=
-
∂
v
c
∂
z
+
γ
x
φ
y
z
=
∂
u
c
∂
z
+
γ
y
,where γ
x and γ
y are the shear strains of y-direction and x-direction due to deflection, respectively.When the curvature is great,φ
x also varies with x and thus it should be written as(2)
φ
x
x
,
z
=
-
∂
v
c
∂
z
1
+
γ
x
.It follows from∂
z
1
=
(
R
0
-
x
/
R
0
)
∂
z that(3)
φ
x
x
,
z
≈
R
0
R
0
-
x
φ
x
z
=
1
+
κ
x
φ
x
z
,which reflects the influence of initial curvature onφ
x and where κ is the curvature of radial direction satisfying κ
=
1
/
R
0.And then, neglecting distortion, we can express the displacement at generic point(
x
,
y
,
z
) within curved box girder in terms of such seven basic displacements, say, u
c
(
z
), v
c
(
z
), w
c
(
z
), φ
x, φ
y, θ
(
z
), and β
(
z
) as (4)
u
x
,
y
,
z
=
u
c
z
-
y
-
y
M
θ
z
v
x
,
y
,
z
=
v
c
z
+
x
θ
z
w
x
,
y
,
z
=
w
c
z
-
φ
y
z
+
κ
w
c
z
x
-
e
+
y
-
y
S
1
+
κ
x
φ
x
z
-
ω
β
′
z
-
κ
φ
x
z
,where e is the eccentricity of flexural center with respect to centroid in the x-direction; ω is the principle sectorial coordinate of the cross section box girder; y
S and y
M are the coordinates of centroid and shear center in y-direction, respectively.


## 3.2. The Relation between the Longitudinal Normal Strain inz-Direction at Generic Point within the Box Girder and Seven Basic Displacements

LettingR
0 and R
0
-
x to be the longitudinal radius of curvature at the centroid (i.e., x
=
0) and generic point of the cross section of box girder, respectively, we can obtain(5)
∂
z
1
=
R
0
-
x
∂
z
R
0
,where ∂
z
1 and ∂
z are the partial derivative of z
1 and z.The abscissa isx, respectively. And herein z
1-axis is parallel to z-axis. Therefore, the longitudinal normal strain in z-direction at generic point within curved box girder considering the effect of initial curvature can be written as [7](6)
ε
z
=
R
0
R
0
-
x
·
∂
w
∂
z
-
u
R
0
-
x
=
1
+
κ
x
∂
w
∂
z
-
κ
u
.Substituting (4) into (6) yields(7)
ε
z
=
1
+
κ
x
w
C
′
-
φ
y
+
κ
w
c
x
-
e
′
+
1
+
κ
x
y
-
y
S
φ
x
′
-
ω
β
′
-
κ
φ
x
′
-
κ
u
c
-
y
-
y
M
θ
.The principle sectorial coordinateω in (7) should be discussed for such cases as top deck, bottom deck, right web, left web, right flange, and left flange, respectively. Due to the effect of variable curvature and variable height, h, κ, y
S, y
M, e, and ω are all functions of z.


## 3.3. The Relations between the Shear Strain at Generic Point within Curved Box Girder and Seven Basic Displacements

Considering any infinitesimal curved element ofd
z by d
s in top flange, bottom flange, or web, we can obtain the shear strain (
γ
) between the plane s
z with curvature taken into account due to a deformation w in the z-direction and u
s:(8)
γ
=
∂
w
∂
s
+
1
+
κ
x
∂
u
s
∂
z
,where u
s (see Figure 2) is the displacement measured along the center line of the closed cross section denoted by s. For the cases of top flange, right web, bottom flange, and left web, u
s is given as (9)
top  flange:
u
s
=
u
x
,
t
1
2
,
z
right  web:
u
s
=
-
u
x
,
y
,
z
cos
⁡
α
+
v
x
,
y
,
z
sin
⁡
α
bottom  flange:
u
s
=
-
u
x
,
h
+
t
1
2
,
z
left  web:
u
s
=
-
u
x
,
y
,
z
cos
⁡
α
-
v
x
,
y
,
z
sin
⁡
α
,where α is the inclination of web.Figure 2
u
s of curved cross section.Substituting (9) into (8) yields the relations of shear stain of top flange, right web, bottom flange, and left web with the displacements, respectively.


## 4. Dynamic Finite Element of Curved Box Girder

### 4.1. The Model of Curved Girder Element and the Corresponding Shape Functions

Set up a 1D generalized local curvilinear coordinate axisζ within the element, in which ζ-direction is the same as z-direction of curvilinear coordinate in global coordinate system and ζ satisfies ζ
=
0 at the middle node and ζ
=
±
1 at two end nodes. And then the element shape functions can be expressed as (10)
N
1
ζ
=
-
1
2
ζ
1
-
ζ
N
2
ζ
=
1
-
ζ
2
N
3
ζ
=
1
2
ζ
1
+
ζ
.The value ofz at generic point within the element in global coordinate system can be obtained by interpolation with shape functions. Meanwhile, it is assumed that κ and h are variable continuously within the element and thus we can obtain κ and h at generic point within element via interpolation similar to that applied to z. Together, they can be written as(11)
z
ζ
=
∑
i
=
1
3
z
i
N
i
ζ
κ
ζ
=
∑
i
=
1
3
κ
i
N
i
ζ
h
ζ
=
∑
i
=
1
3
h
i
N
i
ζ
.Furthermore, we can get Jacobian coefficient:(12)
J
=
d
z
d
ζ
=
∑
i
=
1
3
z
i
d
N
i
ζ
d
ζ
.So far we have established a curved finite element of three nodes that can reflect the variation of curvature within element and variation of height of cross section.


### 4.2. Strain-Displacement Matrix

The nodal displacement of curved beam element considering tension, compression, bending, torsion, and warping can be written as a vector comprised of seven basic displacements:(13)
δ
i
e
=
δ
1
T
δ
2
T
δ
3
T
T
δ
i
=
u
c
i
v
c
i
w
c
i
ϕ
x
i
ϕ
y
i
θ
i
β
i
′
T
i
=
1,2
,
3
.The displacement at generic pointζ can be obtained by interpolation with shape functions in (10) as(14)
δ
=
N
7
×
21
δ
e
,where δ
=
u
c
v
c
w
c
ϕ
x
ϕ
y
θ
β
′
T, N
7
×
21
=
N
1
N
2
N
3, N
i
=
N
i
I
7
×
7
i
=
1,2
,
3, and I is a 7 by 7 identity.According to the strain-displacement relation, we can obtain the normal strain and shear strain as(15)
ε
=
ε
z
γ
=
L
1
L
2
2
×
7
N
7
×
21
δ
e
=
B
2
×
21
δ
ein whichB
2
×
21
=
L
1
L
2
2
×
7
N
7
×
21, and L
1, L
2 can be obtained according to (7) and (8).


### 4.3. Stiffness Matrix, Mass Matrix, and Damping Matrix

#### 4.3.1. Element Stiffness Matrix

According to minimum potential energy principle, the element stiffness matrix can be induced as(16)
K
e
=
∫
-
1
+
1
∬
1
-
κ
x
J
B
T
D
B
d
x
d
y
d
ζ
.Substituting[
B
] in (15) and the stress-displacement matrix [
D
] into (16) will yield its explicit expression and herein it is skipped.


#### 4.3.2. Element Mass Matrix

In this section, we will drive the consistent mass matrix of curved box girder element of three nodes from the kinetic energy expression of element.Taking the first-order derivative of (4) with respect to time t gives the velocity of generic point within the cross section of curved girder element, which is in terms of the derivative of the basic displacements with respect to time t as follows: (17)
u
˙
x
,
y
,
z
=
u
˙
c
z
-
y
-
y
M
θ
˙
z
v
˙
x
,
y
,
z
=
v
˙
c
z
+
x
θ
˙
z
w
˙
x
,
y
,
z
=
w
˙
c
z
-
ϕ
˙
y
z
+
κ
w
˙
c
z
x
-
e
+
y
-
y
S
1
+
κ
x
ϕ
˙
x
z
-
ω
β
˙
′
z
-
κ
ϕ
˙
x
z
.The derivative of the basic displacements within any cross section of element about timet can be expressed by interpolation with shape functions and the derivative of nodal basic displacements within curved girder element about time t as (18)
δ
˙
=
N
7
×
21
δ
˙
e
.Rearranging (18) using (17) gives(19)
u
˙
v
˙
w
˙
=
T
δ
˙
=
T
N
δ
˙
e
,where matrix T is the coefficient matrix derived from (11) and its explicit expression is skipped here.Since the nodal velocities of element have been obtained in (19), the kinetic energy expression in terms of the nodal velocities can be expressed as (20)
1
2
∭
u
˙
v
˙
w
˙
ρ
u
˙
v
˙
w
˙
d
V
=
1
2
∭
δ
˙
e
T
N
T
T
T
ρ
T
N
δ
˙
e
d
x
d
y
R
0
-
x
R
0
d
z
=
1
2
∭
1
-
κ
x
δ
˙
e
T
N
T
T
T
ρ
T
N
δ
˙
e
d
x
d
y
d
z
=
1
2
δ
˙
e
T
∫
-
1
1
∬
1
-
κ
x
J
N
T
T
T
ρ
T
N
d
x
d
y
d
ζ
δ
˙
e
,where ρ is the density of material.From (20) the mass matrix follows as(21)
M
e
=
∫
-
1
1
∬
1
-
κ
x
J
N
T
T
T
ρ
T
N
d
x
d
y
d
ζ
.


#### 4.3.3. Damping Matrix

Generally speaking, the damping of structures cannot be determined exactly. In reality it is obtained approximately by measuring the energy dissipation of structural system during vibration. In doing so, all energy-dissipating mechanisms are taken into account and the global damping matrix is derived straightforward. Unlike the stiffness matrix and mass matrix of structure, it is not necessary to construct the global damping matrix from the element damping matrix by assembling technique and thus no damping matrix of element is needed to be derived. In this paper Rayleigh damping is adopted. It is a linear combination of mass matrix and stiffness matrix as follows:(22)
C
=
α
M
+
β
K
.According to the orthogonality relations between modes with distinct frequencies about mass matrix or stiffness matrix, it follows that Rayleigh damping also satisfies orthogonality relation.


### 4.4. The Solution to the Differential Equation of Motion

For the damping structural system, the differential equation of motion can be written as(23)
M
δ
¨
t
+
C
δ
˙
t
+
K
δ
t
=
F
.Herein mode superposition method is applied to solve (23) for δ
t. To do so, the natural frequencies and modes of structures should be obtained first by solving the frequency equation shown in (24)
-
ω
2
M
+
K
=
0
.A finite element procedure is developed utilizing MATLAB to finish such a task easily. Calling the functionΩ
,
Φ
=
e
i
g
s
K
,
M to solve (20) will yield both natural frequencies and modes. What is more, the modes have been normalized as(25)
Ω
=
diag
⁡
ω
1
,
ω
2
,
…
,
ω
n
Φ
=
φ
1
,
φ
2
,
…
,
φ
n
.Defining the modal expansion of displacement vectorδ
t as δ
t
=
Φ
η
t, and then substituting it into (23), and furthermore multiplying both sides by Φ
T give(26)
Φ
T
M
Φ
η
¨
+
Φ
T
C
Φ
η
˙
+
Φ
T
K
Φ
η
=
Φ
T
F
.It follows from the orthogonality relations of modes in junction with (22) that(27)
φ
i
T
M
φ
j
=
1
i
=
j
0
i
≠
j
φ
i
T
K
φ
j
=
ω
i
2
i
=
j
0
i
≠
j
φ
i
T
C
φ
j
=
α
+
β
ω
i
2
=
2
ξ
i
ω
i
i
=
j
0
i
≠
j
.Substituting them into (26) results in (28)
η
¨
i
+
2
ξ
i
ω
i
η
˙
i
+
ω
i
2
η
i
=
f
i
i
=
1,2
,
3
⋯
n
.Thus, the set ofn coupled differential equations (see (26)) in modal displacements δ
i
t
—
i
=
1,2
,
…
,
n has been transformed to the set of n uncoupled equations (see (28)) in modal coordinates η
i
—
i
=
1,2
,
…
,
n. There are n such equations as (28) and one for each mode. Solving them for η
i
—
i
=
1,2
,
…
,
n and then summing up all modal contributions lead to the total dynamic response:(29)
δ
t
=
∑
i
=
1
n
φ
i
η
i
t
.In reality, the contribution of high order modes to the dynamic response of structure is comparatively negligible and thus truncated. Only the contribution of first several order modes remains, and the rest will not be considered.


## 4.1. The Model of Curved Girder Element and the Corresponding Shape Functions

Set up a 1D generalized local curvilinear coordinate axisζ within the element, in which ζ-direction is the same as z-direction of curvilinear coordinate in global coordinate system and ζ satisfies ζ
=
0 at the middle node and ζ
=
±
1 at two end nodes. And then the element shape functions can be expressed as (10)
N
1
ζ
=
-
1
2
ζ
1
-
ζ
N
2
ζ
=
1
-
ζ
2
N
3
ζ
=
1
2
ζ
1
+
ζ
.The value ofz at generic point within the element in global coordinate system can be obtained by interpolation with shape functions. Meanwhile, it is assumed that κ and h are variable continuously within the element and thus we can obtain κ and h at generic point within element via interpolation similar to that applied to z. Together, they can be written as(11)
z
ζ
=
∑
i
=
1
3
z
i
N
i
ζ
κ
ζ
=
∑
i
=
1
3
κ
i
N
i
ζ
h
ζ
=
∑
i
=
1
3
h
i
N
i
ζ
.Furthermore, we can get Jacobian coefficient:(12)
J
=
d
z
d
ζ
=
∑
i
=
1
3
z
i
d
N
i
ζ
d
ζ
.So far we have established a curved finite element of three nodes that can reflect the variation of curvature within element and variation of height of cross section.


## 4.2. Strain-Displacement Matrix

The nodal displacement of curved beam element considering tension, compression, bending, torsion, and warping can be written as a vector comprised of seven basic displacements:(13)
δ
i
e
=
δ
1
T
δ
2
T
δ
3
T
T
δ
i
=
u
c
i
v
c
i
w
c
i
ϕ
x
i
ϕ
y
i
θ
i
β
i
′
T
i
=
1,2
,
3
.The displacement at generic pointζ can be obtained by interpolation with shape functions in (10) as(14)
δ
=
N
7
×
21
δ
e
,where δ
=
u
c
v
c
w
c
ϕ
x
ϕ
y
θ
β
′
T, N
7
×
21
=
N
1
N
2
N
3, N
i
=
N
i
I
7
×
7
i
=
1,2
,
3, and I is a 7 by 7 identity.According to the strain-displacement relation, we can obtain the normal strain and shear strain as(15)
ε
=
ε
z
γ
=
L
1
L
2
2
×
7
N
7
×
21
δ
e
=
B
2
×
21
δ
ein whichB
2
×
21
=
L
1
L
2
2
×
7
N
7
×
21, and L
1, L
2 can be obtained according to (7) and (8).


## 4.3. Stiffness Matrix, Mass Matrix, and Damping Matrix

### 4.3.1. Element Stiffness Matrix

According to minimum potential energy principle, the element stiffness matrix can be induced as(16)
K
e
=
∫
-
1
+
1
∬
1
-
κ
x
J
B
T
D
B
d
x
d
y
d
ζ
.Substituting[
B
] in (15) and the stress-displacement matrix [
D
] into (16) will yield its explicit expression and herein it is skipped.


### 4.3.2. Element Mass Matrix

In this section, we will drive the consistent mass matrix of curved box girder element of three nodes from the kinetic energy expression of element.Taking the first-order derivative of (4) with respect to time t gives the velocity of generic point within the cross section of curved girder element, which is in terms of the derivative of the basic displacements with respect to time t as follows: (17)
u
˙
x
,
y
,
z
=
u
˙
c
z
-
y
-
y
M
θ
˙
z
v
˙
x
,
y
,
z
=
v
˙
c
z
+
x
θ
˙
z
w
˙
x
,
y
,
z
=
w
˙
c
z
-
ϕ
˙
y
z
+
κ
w
˙
c
z
x
-
e
+
y
-
y
S
1
+
κ
x
ϕ
˙
x
z
-
ω
β
˙
′
z
-
κ
ϕ
˙
x
z
.The derivative of the basic displacements within any cross section of element about timet can be expressed by interpolation with shape functions and the derivative of nodal basic displacements within curved girder element about time t as (18)
δ
˙
=
N
7
×
21
δ
˙
e
.Rearranging (18) using (17) gives(19)
u
˙
v
˙
w
˙
=
T
δ
˙
=
T
N
δ
˙
e
,where matrix T is the coefficient matrix derived from (11) and its explicit expression is skipped here.Since the nodal velocities of element have been obtained in (19), the kinetic energy expression in terms of the nodal velocities can be expressed as (20)
1
2
∭
u
˙
v
˙
w
˙
ρ
u
˙
v
˙
w
˙
d
V
=
1
2
∭
δ
˙
e
T
N
T
T
T
ρ
T
N
δ
˙
e
d
x
d
y
R
0
-
x
R
0
d
z
=
1
2
∭
1
-
κ
x
δ
˙
e
T
N
T
T
T
ρ
T
N
δ
˙
e
d
x
d
y
d
z
=
1
2
δ
˙
e
T
∫
-
1
1
∬
1
-
κ
x
J
N
T
T
T
ρ
T
N
d
x
d
y
d
ζ
δ
˙
e
,where ρ is the density of material.From (20) the mass matrix follows as(21)
M
e
=
∫
-
1
1
∬
1
-
κ
x
J
N
T
T
T
ρ
T
N
d
x
d
y
d
ζ
.


### 4.3.3. Damping Matrix

Generally speaking, the damping of structures cannot be determined exactly. In reality it is obtained approximately by measuring the energy dissipation of structural system during vibration. In doing so, all energy-dissipating mechanisms are taken into account and the global damping matrix is derived straightforward. Unlike the stiffness matrix and mass matrix of structure, it is not necessary to construct the global damping matrix from the element damping matrix by assembling technique and thus no damping matrix of element is needed to be derived. In this paper Rayleigh damping is adopted. It is a linear combination of mass matrix and stiffness matrix as follows:(22)
C
=
α
M
+
β
K
.According to the orthogonality relations between modes with distinct frequencies about mass matrix or stiffness matrix, it follows that Rayleigh damping also satisfies orthogonality relation.


## 4.3.1. Element Stiffness Matrix

According to minimum potential energy principle, the element stiffness matrix can be induced as(16)
K
e
=
∫
-
1
+
1
∬
1
-
κ
x
J
B
T
D
B
d
x
d
y
d
ζ
.Substituting[
B
] in (15) and the stress-displacement matrix [
D
] into (16) will yield its explicit expression and herein it is skipped.


## 4.3.2. Element Mass Matrix

In this section, we will drive the consistent mass matrix of curved box girder element of three nodes from the kinetic energy expression of element.Taking the first-order derivative of (4) with respect to time t gives the velocity of generic point within the cross section of curved girder element, which is in terms of the derivative of the basic displacements with respect to time t as follows: (17)
u
˙
x
,
y
,
z
=
u
˙
c
z
-
y
-
y
M
θ
˙
z
v
˙
x
,
y
,
z
=
v
˙
c
z
+
x
θ
˙
z
w
˙
x
,
y
,
z
=
w
˙
c
z
-
ϕ
˙
y
z
+
κ
w
˙
c
z
x
-
e
+
y
-
y
S
1
+
κ
x
ϕ
˙
x
z
-
ω
β
˙
′
z
-
κ
ϕ
˙
x
z
.The derivative of the basic displacements within any cross section of element about timet can be expressed by interpolation with shape functions and the derivative of nodal basic displacements within curved girder element about time t as (18)
δ
˙
=
N
7
×
21
δ
˙
e
.Rearranging (18) using (17) gives(19)
u
˙
v
˙
w
˙
=
T
δ
˙
=
T
N
δ
˙
e
,where matrix T is the coefficient matrix derived from (11) and its explicit expression is skipped here.Since the nodal velocities of element have been obtained in (19), the kinetic energy expression in terms of the nodal velocities can be expressed as (20)
1
2
∭
u
˙
v
˙
w
˙
ρ
u
˙
v
˙
w
˙
d
V
=
1
2
∭
δ
˙
e
T
N
T
T
T
ρ
T
N
δ
˙
e
d
x
d
y
R
0
-
x
R
0
d
z
=
1
2
∭
1
-
κ
x
δ
˙
e
T
N
T
T
T
ρ
T
N
δ
˙
e
d
x
d
y
d
z
=
1
2
δ
˙
e
T
∫
-
1
1
∬
1
-
κ
x
J
N
T
T
T
ρ
T
N
d
x
d
y
d
ζ
δ
˙
e
,where ρ is the density of material.From (20) the mass matrix follows as(21)
M
e
=
∫
-
1
1
∬
1
-
κ
x
J
N
T
T
T
ρ
T
N
d
x
d
y
d
ζ
.


## 4.3.3. Damping Matrix

Generally speaking, the damping of structures cannot be determined exactly. In reality it is obtained approximately by measuring the energy dissipation of structural system during vibration. In doing so, all energy-dissipating mechanisms are taken into account and the global damping matrix is derived straightforward. Unlike the stiffness matrix and mass matrix of structure, it is not necessary to construct the global damping matrix from the element damping matrix by assembling technique and thus no damping matrix of element is needed to be derived. In this paper Rayleigh damping is adopted. It is a linear combination of mass matrix and stiffness matrix as follows:(22)
C
=
α
M
+
β
K
.According to the orthogonality relations between modes with distinct frequencies about mass matrix or stiffness matrix, it follows that Rayleigh damping also satisfies orthogonality relation.


## 4.4. The Solution to the Differential Equation of Motion

For the damping structural system, the differential equation of motion can be written as(23)
M
δ
¨
t
+
C
δ
˙
t
+
K
δ
t
=
F
.Herein mode superposition method is applied to solve (23) for δ
t. To do so, the natural frequencies and modes of structures should be obtained first by solving the frequency equation shown in (24)
-
ω
2
M
+
K
=
0
.A finite element procedure is developed utilizing MATLAB to finish such a task easily. Calling the functionΩ
,
Φ
=
e
i
g
s
K
,
M to solve (20) will yield both natural frequencies and modes. What is more, the modes have been normalized as(25)
Ω
=
diag
⁡
ω
1
,
ω
2
,
…
,
ω
n
Φ
=
φ
1
,
φ
2
,
…
,
φ
n
.Defining the modal expansion of displacement vectorδ
t as δ
t
=
Φ
η
t, and then substituting it into (23), and furthermore multiplying both sides by Φ
T give(26)
Φ
T
M
Φ
η
¨
+
Φ
T
C
Φ
η
˙
+
Φ
T
K
Φ
η
=
Φ
T
F
.It follows from the orthogonality relations of modes in junction with (22) that(27)
φ
i
T
M
φ
j
=
1
i
=
j
0
i
≠
j
φ
i
T
K
φ
j
=
ω
i
2
i
=
j
0
i
≠
j
φ
i
T
C
φ
j
=
α
+
β
ω
i
2
=
2
ξ
i
ω
i
i
=
j
0
i
≠
j
.Substituting them into (26) results in (28)
η
¨
i
+
2
ξ
i
ω
i
η
˙
i
+
ω
i
2
η
i
=
f
i
i
=
1,2
,
3
⋯
n
.Thus, the set ofn coupled differential equations (see (26)) in modal displacements δ
i
t
—
i
=
1,2
,
…
,
n has been transformed to the set of n uncoupled equations (see (28)) in modal coordinates η
i
—
i
=
1,2
,
…
,
n. There are n such equations as (28) and one for each mode. Solving them for η
i
—
i
=
1,2
,
…
,
n and then summing up all modal contributions lead to the total dynamic response:(29)
δ
t
=
∑
i
=
1
n
φ
i
η
i
t
.In reality, the contribution of high order modes to the dynamic response of structure is comparatively negligible and thus truncated. Only the contribution of first several order modes remains, and the rest will not be considered.


## 5. Numerical Examples

In order to verify the validity and accuracy of the element proposed above, two case studies are performed in this section and it is compared against another two conventional elements.Example 1 (a straight box girder of uniform cross section subjected to a vertical cosine load at midspan).
The project to be studied is a straight box girder of uniform cross section fixed at two ends with span of 33 m subjected to a vertical cosine loadF
=
100
cos
⁡
(
2
π
t
) kN at midspan. The dimensions of cross section are shown in Figure 3. Other given parameters are Young’s modulus E
=
3.25
×
1
0
4 MPa, Poisson’s ratio υ
=
0.1667, and density of material ρ
=
2500 kg/m3. Such a straight box girder can be regarded as a special case of curved box girder by just letting the radius of curvature of the corresponding curved box girder approach infinity. Thus, its dynamic response can be computed in the same way as that applied to curved girder. Three-dimensional beam element BEAM188, a conventional element used in ANSYS software, is also applied herein to analyze the same project such that the validity and accuracy of the present element can be tested by comparison.
The comparison of the frequencies and modes computed by BEAM188 and the element proposed in this paper is listed in Table1. The time-histories of vertical displacement at midspan due to cosine load are shown in Figure 4.
Table1 reveals that the frequencies obtained by the element proposed in this paper agree excellently with those obtained by BEAM188. No error is greater than 8%. The latter element is a conventional element that has been proved effective for straight box girder analysis. As for modes, they are identical to each other. Figure 4 demonstrates that the tendencies of dynamic responses caused by the cosine load at midspan are almost identical for two elements; that is, at the initial phase, they both fluctuates due to transient response, and then transient response decays rapidly and vanishes due to damping, and gradually only steady-state vibration remains.
The comparison above proves that the element proposed in this paper can achieve satisfactory accuracy for the case of straight box girder.Table 1
Comparison of the frequencies and modes.
Mode number
BEAM188/Hz
The present element/Hz
Error/%
The mode that predominates in dynamic response
1
9.372
9.7099
3.61
The first-order mode of vertical bending
2
18.936
18.044
-
4.71
The first-order mode of lateral bending
3
21.717
20.387
-
6.12
The first-order mode of torsion
4
22.177
23.458
5.78
The second-order mode of vertical bending
5
37.845
37.803
-
0.11
The second-order mode of torsion
6
40.317
40.458
0.35
The second-order mode of lateral bending
7
45.197
41.897
-
7.30
The third-order mode of vertical bending
8
54.686
54.630
-
0.10
The first-order mode of axial deformationFigure 3
Dimensions of box girder (unit: m).Figure 4
The time-histories of vertical displacement at midspan.Example 2 (a curved box girder of uniform cross section subjected to a vertical cosine load at midspan).
All given conditions are the same as those of Example1 except that the radius of curvature of center line is 100 m instead of infinity.
Solid element SOLID45 and BEAM188 are also adopted herein to compare with the element proposed in this paper. For curved girder analysis using BEAM188, the curve is approximated by the polygonal lines.
As illustrated in Table2, the results obtained by the element proposed in this paper agree well with those obtained by SOILD45 for some relatively important low order frequencies of vertical bending, especially the first-order frequency. Overall, the errors between the present element and SOILD45 are smaller than those between BEAM188 and SOILD45. And the former are relatively uniform for every frequency and no any mode obviously predominates.
Figure5 shows that the dynamic response obtained by the element in this paper agrees better with that obtained by SOILD45 than that obtained by BEAM188.
The numerical examples above confirm that, for dynamic response analysis, the curved beam element of three nodes developed in this paper can achieve higher accuracy than BEAM188. Compared with SOLID45, it can reduce the complexity of preprocessing and postprocessing and meanwhile show more intuitive results.Table 2
Comparison of the frequencies and modes.
Mode number
①/Hz
②/Hz
③/Hz
Error 1/%
Error 2/%
The mode that dominates dynamic response
1
9.67
9.34
9.69
-
3.39
0.23
The first-order mode of vertical bending
2
19.47
19.42
19.28
-
0.27
-
0.98
The first-order mode of lateral bending
3
21.66
21.89
21.43
1.10
-
1.06
The first-order mode of torsion
4
22.98
22.14
23.42
-
3.66
1.92
The second-order mode of vertical bending
5
36.57
37.81
35.71
3.40
-
2.35
The second-order mode of torsion
6
38.67
40.41
40.11
4.49
3.72
The second-order mode of lateral bending
7
43.74
44.98
42.83
2.83
-
2.08
The third-order mode of vertical bending
8
51.39
55.15
54.91
7.30
6.85
The first-order mode of axial deformation
Notes: ① for SOLID45, ② for BEAM188, and ③ for the present element, error 1 for the difference between the frequencies obtained by BEAMA188 and those obtained by SOILD45, error 2 for the difference between the frequencies obtained by the present element and those obtained by SOILD45.Figure 5
The time-histories of vertical displacement at midspan.


## 6. Conclusions

A finite element of three nodes and twenty-one degrees of freedom for dynamic analysis of curved box girder with such deformations as tension, compression, bending, torsion, and warping taken into account was developed. The relation of normal stress and shear strain at generic point within cross section with the basic displacements of cross section caused by bending and torsion was first established, and then from which the element stiffness matrix was derived. The relation between the velocities at generic point within the element and the velocities of nodal velocity of cross section was first obtained, and then the expression of kinetic energy was yielded, and furthermore the element mass matrix was deduced from such an expression. Rayleigh damping matrix was obtained easily since it is a linear combination of stiffness matrix and mass matrix and the latter two have been derived. Solving the frequency equations of undamped free vibration yielded the natural frequencies of structure. And then solving the differential equations of motion with mode superposition method gave the dynamic responses of structure. A MATLAB procedure was developed to implement such a finite element method.The finite element purposed herein fully takes account of deformation characteristics of curved thin-walled box girder including tension, compression, torsion, and warping. Meanwhile, it can accommodate well to the vibration of curvature within element and the vibration of height of cross section. It can greatly reduce the complexity of preprocessing and postprocessing. Compared with the solid element and plate element, it is more effective and intuitive, especially for dynamic analysis. It can achieve higher accuracy in comparison with the straight beam element.



---

*Source: 8460130-2016-05-11.xml*