# Corotational modeling and NURBS-based kinematic constraint implementation in
three-dimensional vehicle-track-structure interaction analysis

Maria Fedorova111Current address: LARSA, Inc., 68 S Service Rd, Ste 100,
Melville, 11747, NY, USA<EMAIL_ADDRESS>M.V. Sivaselvan<EMAIL_ADDRESS>
###### Abstract

An algorithm for three-dimensional dynamic vehicle-track-structure interaction
(VTSI) analysis is described in this paper. The algorithm is described in
terms of bridges and high-speed trains, but more generally applies to
multibody systems coupled to deformable structures by time-varying kinematic
constraints. Coupling is accomplished by a kinematic constraint/Lagrange
multiplier approach, resulting in a system of index-3 Differential Algebraic
Equations (DAE). Three main new concepts are developed. (i) A corotational
approach is used to represent the vehicle (train) dynamics. Reference
coordinate frames are fitted to the undeformed geometry of the bridge. While
the displacements of the train can be large, deformations are taken to be
small within these frames, resulting in linear (time-varying) rather than
nonlinear dynamics. (ii) If conventional finite elements are used to
discretize the track, the curvature is discontinuous across elements (and
possibly rotation, too, for curved tracks). This results in spurious numerical
oscillations in computed contact forces and accelerations, quantities of key
interest in VTSI. A NURBS-based discretization is employed for the track to
mitigate such oscillations. (iii) The higher order continuity due to using
NURBS allows for alternative techniques for solving the VTSI system. First,
enforcing constraints at the acceleration level reduces an index-3 DAE to an
index-1 system that can be solved without numerical dissipation. Second, a
constraint projection method is proposed to solve an index-3 DAE system
without numerical dissipation by correcting wheel velocities and
accelerations. Moreover, the modularity of the presented algorithm, resulting
from a kinematic constraint/Lagrange multiplier formulation, enables ready
integration of this VTSI approach in existing structural analysis and finite
element software. A simplified vehicle-bridge model is first introduced that
exemplifies numerical challenges. This model is used throughout the paper to
motivate various algorithmic choices. Finally, a numerical example consisting
of realistic vehicle and bridge models is presented.

###### keywords:

Vehicle-track-structure interaction (VTSI) , High-speed railway bridge ,
Differential algebraic equations (DAE) , Lagrange multipliers , Corotational ,
Non-uniform rational B-splines (NURBS)

[inst1]organization=Department of Civil, Structural, and Environmental
Engineering, University at Buffalo,addressline=Ketter Hall, city=Buffalo,
postcode=14260, state=NY, country=USA

##  1 Introduction

Traditional static analysis with moving loads is not sufficient for many
bridges supporting high-speed trains. Typically used codes of practice [1, 2]
require dynamic analysis for such structures. It can be performed either with
a series of moving loads representing the dynamic effects of multiple trains
and axles or by solving a coupled system of vehicle, track, and structure
equations of motion. The latter approach constitutes vehicle-track-structure
interaction (VTSI) analysis. According to the California High-Speed Train
Project documentation [1], VTSI analysis is required for some structures that
do not meet the requirements of a simplified analysis. [3]

There are numerous ways to model dynamic interaction between trains and track-
structure subsystems (see, for example, References [4, 5, 6] for an overview
of available approaches). As was discussed by Fedorova and Sivaselvan [7], one
of the drawbacks of many existing algorithms is the way the governing
equations of motion are solved. The equations for the train and bridge
subsystems are often combined together into a single equation to obtain the
solution [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
25, 26, 27, 28, 29, 30, 31, 32, 33], or the two sets of equations are solved
separately employing iterative procedures [34, 35, 36, 37, 38, 39, 40, 41, 42,
43, 44]. While the former approach makes it difficult to integrate with
existing structural analysis software, the latter requires significant
computational effort. Another possible approach is to solve the system of
equations using a complementary kinematic constraint condition [45, 46, 7].
This approach facilitates integration into existing software with minimum
intervention, and, moreover, modeling of contact separation between the wheels
and the bridge becomes more straightforward. Furthermore, wheel-rail contact
forces are direct outputs of the computation.

This paper builds on an algorithm for two-dimensional vehicle-track-structure
interaction (VTSI) analysis [7]. There are three main new features of the
approach presented in this paper: 1) a corotational approach, wherein train
dynamics is described with reference to coordinate frames fitted to the bridge
geometry; deformations within these frames can be taken to be small, resulting
in linear (time-varying) dynamics; 2) NURBS-based discretization of the
vehicle path to provide sufficient inter-element continuity, and prevent
spurious oscillations in contact forces and accelerations, quantities of
principal interest in VTSI; 3) two alternative ways of solving index-3 system
of differential-algebraic equations (DAE) that are possible due to higher
continuity provided by NURBS: enforcement of constraints at acceleration
level, which reduced the system to a index-1 DAE, and constraint projection
method, which allows for solving an index-3 system without numerical
dissipation by correcting wheel velocities and accelerations. Moreover, while
the algorithm is applied specifically to VTSI analysis in the current paper,
it can generally be used in multibody systems coupled to deformable structures
by time-varying kinematic constraints.

This paper is organized as follows. Section 2 introduces the corotational
approach used to describe vehicle kinematics. A numerical example is then
introduced (Section 3) and used as a running thread throughout the paper to
illustrate the essential features of the algorithm and describe numerical
issues and resolutions. These numerical issues are described in Section 4,
associated with solving VTSI equations as an index-3 DAE and using finite
elements to model the track. To resolve the numerical issues, NURBS-based
discretization of the track is employed in Section 5. In Section 6, we
describe three ways of solving VTSI equations: (i) first, a
Generalized-$\alpha$ scheme is used to obtain desired numerical damping and
mitigate the spurious oscillations; (ii) second, the advantage of using NURBS
for modeling the track is demonstrated by explicitly enforcing constraints at
acceleration level and avoiding the use of numerical dissipation; (iii) third,
we discuss a simple constraint projection method that allows solving VTSI as
an index-3 DAE system without numerical dissipation by non-iterative
correction of obtained velocities and accelerations at each time step. The
advantages and disadvantages of the proposed methods in terms of practical
implementation are discussed. While some of the methods are easy to implement,
they may not be well suited for existing structural analysis software. Further
numerical results for a simplified model, as well as sample results for a
realistic model, are included in Section 7. The key findings are then
summarized in Section 8.

##  2 Corotational approach

The corrotational approach comes from continuum mechanics and the finite
element method (see, for example, Felippa et al. [47] for an overview of the
approach and [48] for an example of practical implementation). The key idea of
this approach is to allow the displacements and rotations of a reference
coordinate frame to be large while keeping deformations within the frame
small. In the VTSI context, cars of trains transversing curved paths undergo
large rigid body motion, but deformations, for instance, of the suspension
system, are small. It is, therefore, convenient to express the car dynamics in
local corotated coordinate frames fit to the curved profile of the path. In
this paper, reference Frenet frames are fit to the path of the vehicle (Figure
1). At the start of the analysis, the vehicle’s initial position at the
beginning of the path is described with respect to the initial Frenet frame
(see Figure 1). Later, as the vehicle moves, the deformations within each car
of the train are computed in a Frenet frame corresponding to the position of
the car body’s center of gravity along the curved path of the bridge. Since
the train longitudinal dynamics are not modeled in this study, longitudinal
connections between the cars (couplers) are not modeled either. The
displacements and rotations of the car parts inside the Frenet frame are taken
to be small. Each part is taken to be a rigid body, forming a system of rigid
bodies (see Figure 2 for a simplified vehicle model and Figure 4 for a
realistic train model). As opposed to trajectory coordinates often used in
railroad vehicle formulations (for example, Shabana et al. [5]), the
corotational approach requires deformations within the reference frames to be
small, resulting in linear dynamics, which is appropriate in practical
situations, but simplifies the VTSI formulation.

Figure 1: The corotational approach in VTSI: motion of a body is described
with respect to a moving Frenet frame. Displacements and rotations of the
Frenet frame can be large, but displacements and rotations relative to the
Frenet frame are small. Here, $XYZ$ is a Global frame, $tnb$ is a Frenet
frame, and $xyz$ is a rigid body frame. The apostrophe indicates the frames
moved to a new, corotated position.

The path taken by a train is a curve in 3D. Let
${\bf\hat{x}}\mathrel{\mathop{\mathchar
58\relax}}[0,L]\rightarrow\mathbb{R}^{3}$ be such a curve parametrized by arc
length $s$ (in Section 5, this curve will be represented using NURBS). Using
parameter $\xi=\hat{\xi}(s)$, this curve can be reparameterized as
${\bf\tilde{x}}\mathrel{\mathop{\mathchar
58\relax}}[\xi_{a},\xi_{b}]\rightarrow\mathbb{R}^{3}$, where the end points
$\xi_{a}$ and $\xi_{b}$ correspond to the start and end of the path. A local
Frenet frame $\\{{\bf t},{\bf n},{\bf b}\\}$ can then be defined at each point
${\bf\tilde{x}}(\xi)$ of the curve as [49, 50]

$\textbf{t}=\frac{{\bf\tilde{x}}^{\prime}}{\left\|{\bf\tilde{x}}^{\prime}\right\|},\
\textbf{n}=\textbf{b}\times\textbf{t},\
\textbf{b}=\frac{{\bf\tilde{x}}^{\prime}\times{\bf\tilde{x}}^{\prime\prime}}{\left\|{\bf\tilde{x}}^{\prime}\times{\bf\tilde{x}}^{\prime\prime}\right\|}$
(1)

where a prime denotes a derivative with respect to parameter $\xi$;
$\textbf{t},\textbf{n},$ and b are the tangent, normal, and binormal vectors.
The Frenet frame is attached to the curve and serves both as a reference frame
to describe the train dynamics (Section 3.1) and as a basis for the geometry
and solution space of equations of motion (Section 5) of a beam representing a
simplified bridge/track (Section 3.2) or rails in the realistic model (Section
3.4).

Having introduced the corotational approach in the VTSI context, we next
develop the numerical examples that are used to motivate the algorithmic
choices.

##  3 Numerical example to motivate algorithmic choices

###  3.1 Simplified vehicle model

In developing our VTSI algorithm, a number of choices are encountered along
the way, such as the use of numerical damping to mitigate the spurious
oscillations (Section 6.1), explicit enforcement of the acceleration
constraints (Section 6.2), and the employment of a simple constraint
projection method (Section 6.3). These choices are best explained by means of
an example. Here, a simplified model is introduced for this purpose. The model
is used to illustrate essential features of the three-dimensional VTSI
algorithm and numerical issues associated with solving an index-3 DAE system.
The model consists of two rigid bodies representing a wheel and a car and a
spring and dashpot representing a vehicle suspension, as shown in Figure 2.
The bridge itself is represented in the simplified model by one curved beam,
without explicit rails, sleepers, ballast etc. (Section 3.2). Application of
the algorithm to realistic models including such features is demonstrated in
Section 7.2.

####  3.1.1 Kinematics

The simplified vehicle model has 4 degrees of freedom. Therefore, four
generalized coordinates are required to describe the motion of the model: the
transverse $u^{\text{t}}_{1}$ and vertical $u^{\text{t}}_{2}$ displacements of
the wheel, rolling of the vehicle $u^{\text{t}}_{3}$, and vertical
displacement of the car $u^{\text{t}}_{4}$. The superscript ”t” stands for
”train”. Figure 2 shows the simplified vehicle in the Frenet frame attached to
the undeformed central axis of the bridge.

Figure 2: Bridge curved in plan and a simplified vehicle model (generalized
coordinates of the vehicle model with respect to the Frenet frame are shown in
blue; position vectors of the vehicle parts are shown in red)

The positions of the simplified wheel (${\bf x}^{\text{w}}$) and car (${\bf
x}^{\text{c}}$) relative to the global frame are then

$\displaystyle{\bf x}^{\text{w}}$ $\displaystyle={\bf x}^{\text{F}}+{\bf
R}^{\text{F}}\left({\begin{array}[]{*{20}{c}}0\\\ {u^{\text{t}}_{1}}\\\
{u^{\text{t}}_{2}}\end{array}}\right)={\bf x}^{\text{F}}+{\bf
R}^{\text{F}}{\begin{bmatrix}{}0&0&0&0\\\ 1&0&0&0\\\
0&1&0&0\end{bmatrix}}\left({\begin{array}[]{*{20}{c}}u^{\text{t}}_{1}\\\
u^{\text{t}}_{2}\\\ u^{\text{t}}_{3}\\\
u^{\text{t}}_{4}\end{array}}\right)={\bf x}^{\text{F}}+{\bf R}^{\text{F}}{\bf
T}^{\text{w}}{\bf u}^{\text{t}}$ (9) $\displaystyle{\bf x}^{\text{c}}$
$\displaystyle={\bf x}^{\text{F}}+{\bf
R}^{\text{F}}\left({\begin{array}[]{c}0\\\
{u^{\text{t}}_{1}-u^{\text{t}}_{3}{l_{0}}}\\\
{u^{\text{t}}_{4}+{l_{0}}}\end{array}}\right)={\bf x}^{\text{F}}+{\bf
R}^{\text{F}}{\begin{bmatrix}{}0&0&0&0\\\ 1&0&-l_{0}&0\\\
0&0&0&1\end{bmatrix}}\left({\begin{array}[]{{c}}u^{\text{t}}_{1}\\\
u^{\text{t}}_{2}\\\ u^{\text{t}}_{3}\\\
u^{\text{t}}_{4}\end{array}}\right)+{\bf
R}^{\text{F}}\left({\begin{array}[]{{c}}0\\\ 0\\\ l_{0}\end{array}}\right)$
(20) $\displaystyle={\bf x}^{\text{F}}+{\bf R}^{\text{F}}{\bf
T}^{\text{c}}{\bf u}^{\text{t}}+{\bf
R}^{\text{F}}\left({\begin{array}[]{c}0\\\ 0\\\ l_{0}\end{array}}\right)$ (24)

where ${\left({\bf x}^{\text{F}},{\bf R}^{\text{F}}\right)}$ is the
specification of the configuration of the Frenet frame relative to the global
frame: ${\bf x}^{\text{F}}$ is the position vector of the origin of the frame
and ${\bf R}^{\text{F}}$ is an orientation of the frame.

The velocities of the wheel and the car in the global frame are,
differentiating (9),

$\displaystyle{\bf\dot{x}}^{\text{w}}$
$\displaystyle={\bf\dot{x}}^{\text{F}}+{\bf
R}^{\text{F}}{\boldsymbol{\hat{\omega}}}^{\text{F}}{\bf T}^{\text{w}}{\bf
u}^{\text{t}}+{\bf R}^{\text{F}}{\bf T}^{\text{w}}\dot{\bf u}^{\text{t}}$ (25)
$\displaystyle{\bf\dot{x}}^{\text{c}}$
$\displaystyle={\bf\dot{x}}^{\text{F}}+{\bf
R}^{\text{F}}{\boldsymbol{\hat{\omega}}}^{\text{F}}{\bf T}^{\text{c}}{\bf
u}^{\text{t}}+{\bf R}^{\text{F}}{\bf T}^{\text{c}}\dot{\bf u}^{\text{t}}+{\bf
R}^{\text{F}}{\boldsymbol{\hat{\omega}}}^{\text{F}}\left({\begin{array}[]{c}0\\\
0\\\ l_{0}\end{array}}\right)$

where ${\boldsymbol{\hat{\omega}}}^{\text{F}}={{\bf
R}^{\text{F}}}^{T}{\bf\dot{R}}^{\text{F}}$ is the matrix representing the
angular velocity of the Frenet frame [51]. The matrix
${\boldsymbol{\hat{\omega}}}^{\text{F}}$ is skew-symmetric by definition and
serves to map $\mathbb{R}^{3}\mapsto\mathbb{R}^{3\times 3}$, that is, for the
angular velocity vector
${\boldsymbol{\omega}}^{\text{F}}={\left({\omega}^{\text{F}}_{1}\
{\omega}^{\text{F}}_{2}\ {\omega}^{\text{F}}_{3}\right)}^{T}$, the matrix
${\boldsymbol{\hat{\omega}}}^{\text{F}}$is defined as

${\boldsymbol{\hat{\omega}}}^{\text{F}}=\left[{\begin{array}[]{*{20}{c}}0&-{\omega}^{\text{F}}_{3}&{\omega}^{\text{F}}_{2}\\\
{\omega}^{\text{F}}_{3}&0&-{\omega}^{\text{F}}_{2}\\\
-{\omega}^{\text{F}}_{2}&{\omega}^{\text{F}}_{1}&0\end{array}}\right]$ (26)

Using properties of cross product operation between vectors [51], we then can
write

${\boldsymbol{\omega}}^{\text{F}}\times{\bf{x}}={\boldsymbol{\hat{\omega}}}^{\text{F}}{\bf{x}}$
(27)

Having described vehicle kinematics, we proceed to deriving the equations of
motion.

####  3.1.2 Equations of motion

The motion of the system is described by employing Lagrangian mechanics and
using the Euler-Lagrange equations of motion for a constrained system [51]
(summation over the index $j$ is assumed):

$\frac{d}{{dt}}\left({\frac{{\partial
L}}{{\partial{{\dot{u}}_{i}}}}}\right)-\frac{{\partial
L}}{{\partial{u_{i}}}}+\frac{\partial C_{j}}{\partial u_{i}}\lambda_{j}=0$
(28)

where ${L}$ is a Lagrangian for an unconstrained system, ${u_{i}}$ are
generalized coordinates of the model, $C_{j}(u)=0$ is a constraint (35c), and
$\lambda_{j}$ is the corresponding Lagrange multiplier. The Lagrangian is
constructed from the kinetic energy, ${T}$, and the potential energy, ${V}$,
as follows

$L=T-V$ (29)

By substituting (29) into (28), and taking into account the fact that the
potential energy ${V}$ is not a function of time or velocity, we obtain

$\frac{d}{{dt}}\left({\frac{{\partial
T}}{{\partial{{\dot{u}}_{i}}}}}\right)-\frac{{\partial
T}}{{\partial{u_{i}}}}+\frac{{\partial V}}{{\partial{u_{i}}}}+\frac{\partial
C_{j}}{\partial u_{i}}\lambda_{j}=0\\\ $ (30)

The kinetic and potential energies of the simplified train model are given by

$\begin{gathered}T=\frac{1}{2}m^{\text{w}}{{{\bf{\dot{x}}}}^{{{\text{w}}^{\mathsf{T}}}}}{\bf\dot{x}}^{\text{w}}+\frac{1}{2}m^{\text{c}}{{{\bf{\dot{x}}}}^{{{\text{c}}^{\mathsf{T}}}}}{\bf\dot{x}}^{\text{c}}+\frac{1}{2}\left(I^{\text{w}}+I^{\text{c}}\right){\dot{u}_{3}}^{{\text{t}^{2}}}\\\
V=\frac{1}{2}k^{\text{s}}(u^{\text{t}}_{4}-u^{\text{t}}_{2})^{2}+m^{\text{w}}g\left(x^{\text{w}}_{3}-x^{\text{w}}_{3}\left(0\right)\right)+m^{\text{c}}g\left(x^{\text{c}}_{3}-x^{\text{c}}_{0}\left(0\right)\right)\\\
\end{gathered}$ (31)

where $m^{\text{w}}$, $m^{\text{c}}$, $I^{\text{w}}$, and $I^{\text{c}}$ are
masses and moments of inertia of the simplified wheel and car, and
$k^{\text{s}}$ is the stiffness of the suspension spring.

By substituting (25) into (31), and from the Euler-Lagrange equation (30), we
obtain the equations of motion of the model

${\bf M}^{\text{t}}\ddot{\bf u}^{\text{t}}+{\bf C}^{\text{t}}(t)\dot{\bf
u}^{\text{t}}+{\bf K}^{\text{t}}(t){\bf u}^{\text{t}}+{{\bf
L}^{\text{t}}}^{\mathsf{T}}\boldsymbol{\lambda}(t)={\bf P}^{\text{t}}(t)$ (32)

where the train matrices are stated explicitly in Appendix A. The complete
derivation of Equation (32) can be found in Fedorova [52].

In the numerical example used throughout this paper to motivate algorithmic
choices, the vehicle properties are partially adapted from a numerical example
by Yang et al. [4]. The element representing a car of the simplified vehicle
model has the mass ${m^{\text{c}}=41750\ $\mathrm{k}\mathrm{g}$}$, mass moment
of inertia $I^{\text{c}}_{x}=23.2\times 10^{3}\
$\mathrm{k}\mathrm{g}\cdot\mathrm{m}^{2}$$, and its center of gravity located
at $1.37\ $\mathrm{m}$$ above the wheel. The wheel element has a mass equal to
the mass of four wheelsets ${m^{\text{w}}=1780\times 4=7120\
$\mathrm{k}\mathrm{g}$}$, and mass moment of inertia
$I^{\text{c}}_{w}=1.14\times 10^{3}\
$\mathrm{k}\mathrm{g}\cdot\mathrm{m}^{2}$$. The listed mass moments of inertia
$I^{\text{c}}_{x}$ and $I^{\text{w}}_{x}$ belong to the Shinkansen (SKS)
series 300 vehicle model[4]. The stiffness of the suspension spring is the
equivalent stiffness of all primary and secondary suspension springs of a
train car
${k^{\text{s}}=2/\left(\frac{1}{4k_{\text{prim}}}+\frac{1}{2k_{\text{sec}}}\
\right)=865.6\ $\mathrm{k}\mathrm{N}\mathrm{/}\mathrm{m}$}$. Here, two bogies
per car and two wheelsets per bogie are assumed, with two springs per each
suspension. The initial length of the spring is ${l_{0}=1.37\ $\mathrm{m}$}$
and corresponds to the position of the car body’s center of gravity. The speed
of the train ${v^{\text{t}}}$ is constant and equals to ${360\
$\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{h}$=100\
$\mathrm{m}\mathrm{/}\mathrm{s}$}$, which equals the maximum speed of the real
trains used in Eurocode EN 1991-2:2003 [2].

###  3.2 Simplified bridge model

A simplified bridge model is first used for clarity of presentation; a
realistic model is used in the numerical example in Section 7.2. The
simplified model does not distinguish between the track and the deck (Figure
2). The bridge can be modeled using standard FEM elements or using NURBS
(Section 5). The numerical example used throughout this paper employs a
simplified bridge model that consists of two straight spans, two transition
spans, and one curved span in the middle of the structure (Figure 3a). The
bridge is fixed at the ends and simply-supported in the middle supports.

The material and geometric properties of the bridge are partially adapted from
a straight bridge model by Yang et al. [4]. The length of each span is ${30\
$\mathrm{m}$}$. According to the Technical Memorandum TM 2.1.2 [53] by the
California High-Speed Train Project, the ”exceptional” (most minimal) track
radius for speed ${355\ $\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{h}$}$ should be
${6000\ $\mathrm{m}$}$. At the same time, the ”desirable” radius for such
speed is ${10700\ $\mathrm{m}$}$. Considering the train speed of ${360\
$\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{h}$}$ (Section 3.1), for the simplified
bridge model the radius of ${6000\ $\mathrm{m}$}$ is chosen. Transition spans
are modeled between the straight and curved spans. The radius of the
transition spans is ${0\ $\mathrm{m}$}$ at the beginning of the span and
$6000\ $\mathrm{m}$$ at the end, with linear variation along the length of the
span. Track with the radius of ${11000\ $\mathrm{m}$}$ is used in the
realistic bridge model (Figure 3b, Sections 3.4 and 7.2).

The following material properties are used: Young’s modulus ${E^{b}=28.25\
$\mathrm{G}\mathrm{P}\mathrm{a}$}$, modulus of rigidity ${G^{b}=10^{3}\
$\mathrm{G}\mathrm{P}\mathrm{a}$}$, torsional constant ${J=15.65\
$\mathrm{m}^{4}$}$, cross-section area ${A=7.73\ $\mathrm{m}^{2}$}$, cross-
section second moments of inertia ${I_{y}^{b}=7.84\ $\mathrm{m}^{4}$}$ and
${I_{z}^{b}=74.42\ $\mathrm{m}^{4}$}$, and the mass per unit length (including
mass of the ballast) ${41740\ $\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{m}$}$.
The standard finite element (constructed with Hermitian cubic shape functions)
or NURBS of degrees $p=3$ and $p=5$ are used.

|
---|---

(a) Simplified bridge model (Section 3.2)

(b) Realistic bridge/track model, ballast elements are not shown (Section 3.4)

Figure 3: Paths in simplified and realistic bridge models used in numerical
examples

###  3.3 Realistic vehicle model

While the simplified vehicle model (Figure 2) in Section 3.1 is meant to
illustrate the essential features of the proposed algorithm, a realistic train
model is also presented for practical three-dimensional VTSI analysis. In this
paper, the use of the realistic model is restricted to Section 7.2. Additional
analysis results and the model derivations are available in Fedorova [54].

The vehicle equations of motion describe the behavior of rigid bodies
(representing train cars, bogies, and wheelsets) and include damping and
stiffness terms corresponding to the suspension system that connects these
rigid bodies

$\begin{gathered}\operatorname*{{{{{\mathcal{A}}}}}}\limits_{i=1}^{n}\left[{\bf
M}^{\text{rb},\ i}{\bf{\ddot{q}}}^{i}+{\bf C}^{\text{rb},\
i}(t){\bf{\dot{q}}}^{i}+{\bf K}^{\text{rb},\ i}(t){\bf
q}^{i}\right]+\operatorname*{{{{{\mathcal{A}}}}}}\limits_{j=1}^{l}\left[{\bf
C}^{\text{susp},\ j}\left(\begin{array}[]{@{}c@{}}{\bf\dot{q}}^{j_{1}}\\\
{\bf\dot{q}}^{j_{2}}\end{array}\right)\right]+\\\
\operatorname*{{{{{\mathcal{A}}}}}}\limits_{k=1}^{m}\left[{\bf
K}^{\text{susp},\ k}\left(\begin{array}[]{@{}c@{}}{\bf\dot{q}}^{k_{1}}\\\
{\bf\dot{q}}^{k_{2}}\end{array}\right)\right]+\operatorname*{{{{{\mathcal{A}}}}}}\limits_{p=1}^{q}\left[{\bf
L}^{\text{rw},\ p}\boldsymbol{\lambda}^{\text{rw},\
p}(t)\right]+\operatorname*{{{{{\mathcal{A}}}}}}\limits_{r=1}^{s}\left[{\bf
L}^{\text{t},\ r}\boldsymbol{\lambda}^{\text{t},\
r}(t)\right]=\operatorname*{{{{{\mathcal{A}}}}}}\limits_{i=1}^{n}\left[{\bf
P}^{\text{rb},\ i}(t)\right]\end{gathered}$ (33)

where ${\bf M}^{\text{rb},\ i}$, ${\bf C}^{\text{rb},\ i}$, ${\bf
K}^{\text{rb},\ i}$ are the mass, damping, and stiffness matrices, ${\bf
P}^{\text{rb},\ i}$ is the load vector, and ${\bf q}^{i}$ is the displacement
vector of $i$-th rigid body; ${\bf C}^{\text{susp},\ j}$ is a damping matrix
and ${\bf\dot{q}}^{j_{1}}$ is a velocity vector of the first joint of $j$-th
suspension dashpot; ${\bf K}^{\text{susp},\ k}$ is a stiffness matrix and
${\bf\dot{q}}^{k_{1}}$ is a velocity vector of the first joint of $k$-th
suspension spring; ${\bf L}^{\text{rw},\ p}$ and
$\boldsymbol{\lambda}^{\text{rw},\ p}$ are wheel-rail constraint matrix and
contact forces of $p$-th wheelset; ${\bf L}^{\text{t},\ r}$ and
$\boldsymbol{\lambda}^{\text{t},\ r}$ are kinematic constraint matrix and
corresponding constraint forces of $r$-th rigid body; $n$ is the number of
rigid bodies (train cars, bogies, and wheelsets), $l$ is the number of
dashpots, $m$ is the number of springs, $q$ is the number of wheelsets, $s$ is
the number of rigid bodies subjected to kinematic constraints;
$\operatorname*{{{{{\mathcal{A}}}}}}$ is the assembly operator.

Equation 33 can be written in a more concise form as

$\begin{gathered}{\bf M}^{\text{t}}{\bf{\ddot{q}}}^{\text{t}}+{\bf
C}^{\text{t}}(t){\bf{\dot{q}}}^{\text{t}}+{\bf K}^{\text{t}}(t){\bf
q}^{\text{t}}+{\bf L}^{\text{rw}}\boldsymbol{\lambda}^{\text{rw}}(t)+{\bf
L}^{\text{t}}\boldsymbol{\lambda}^{\text{t}}(t)={\bf
P}^{\text{t}}(t)\end{gathered}$ (34)

where ${\bf L}^{\text{rw}}$ is a kinematic constraint matrix used to constrain
wheels to the rails and $\boldsymbol{\lambda}^{\text{rw}}$ is a vector of
corresponding contact forces; ${\bf L}^{\text{t}}$ is kinematic constraint
matrix used to restrict longitudinal motions of rigid bodies and rotations of
wheelsets with respect to transverse axis, $n$, and
$\boldsymbol{\lambda}^{\text{t}}$ is a vector of corresponding constraint
forces. An example train car with two bogies and two wheelsets per bogie then
has a total of 31 degrees of freedom, as shown in Figure 4.

Figure 4: An example of a realistic train model (one car is shown)

###  3.4 Realistic bridge model

A realistic bridge model represents the deck and rails as distinct components
and accounts for the ballast (see Figure 3b, where ballast elements are not
shown, and Figure 5). When the analysis is implemented in finite element
analysis software, the deck is modeled with standard finite elements, and the
ballast is modeled with springs acting in vertical and transverse directions.
The rails can be modeled either with NURBS elements (Section 5) or with
standard finite elements. The rationale behind modeling the rails with NURBS
is given in Section 4.2. The Kinematic constraints are used to connect the
rails to the ballast. This approach provides the desired modularity since the
track and the bridge can be modeled independently. For a detailed description
of this model, see Fedorova [54].

Figure 5: Schematic of a realistic bridge and track model

The numerical examples presented in this section are used to motivate
algorithmic choices throughout the paper. In the next section, after
presenting an overview of VTSI equations, we describe encountered numerical
challenges using the simplified VTSI model.

##  4 Summary of VTSI governing equations and numerical challenges

###  4.1 Equations of motion

In this section, the simplified vehicle (Section 3.1) and the simplified
bridge (Section 3.2) are considered. If a realistic VTSI model is considered,
the equations of motion (35) are appended with equations of motion of the
track and kinematic constraints between the track and the bridge. Combining
the equations of motion of the train (32) and the bridge with the kinematic
constraint equation, we obtain the coupled equations of motion

$\displaystyle{\bf M}^{\text{t}}\ddot{\bf u}^{\text{t}}+{\bf
C}^{\text{t}}(t)\dot{\bf u}^{\text{t}}+{\bf K}^{\text{t}}(t){\bf
u}^{\text{t}}+{{\bf L}^{\text{t}}}^{{\mathsf{T}}}\boldsymbol{\lambda}$
$\displaystyle={\bf P}^{\text{t}}(t)$ (35a) $\displaystyle{\bf
M}^{\text{b}}\ddot{\bf u}^{\text{b}}+{\bf C}^{\text{b}}\dot{\bf
u}^{\text{b}}+{\bf K}^{\text{b}}{\bf u}^{\text{b}}+{{\bf
L}^{\text{b}}(t)}^{{\mathsf{T}}}\boldsymbol{\lambda}$ $\displaystyle={\bf
P}^{\text{b}}$ (35b) $\displaystyle{\bf L}^{\text{t}}{\bf u}^{\text{t}}+{\bf
L}^{\text{b}}(t){\bf u}^{\text{b}}$ $\displaystyle=\mathbf{0}$ (35c)

where superscripts ”t”, and ”b” stand for train and bridge; ${\bf
M^{\square}}$, ${\bf C^{\square}}$ and ${\bf K^{\square}}$ are the mass,
damping, and stiffness matrices; ${\bf P^{\square}}$ is the external load
vector; ${{\bf u}^{\square}}$, ${\dot{\bf u}^{\square}}$ and ${\ddot{\bf
u}^{\square}}$ are vectors of displacements, velocities, and accelerations;
$\boldsymbol{\lambda}$ is the vector of contact forces between the train
wheels and the bridge. Matrices ${\bf L}^{\text{t}}$ and ${\bf
L}^{\text{b}}(t)$ are kinematic constraint matrices, where ${\bf
L}^{\text{b}}(t)$ is a time-dependent matrix derived based on the current
locations of train wheels on the bridge. Matrices ${\bf C}^{\text{t}}(t)$ and
${\bf K}^{\text{t}}(t)$ are time-varying due to the corotational formulation.
The vector $\boldsymbol{\lambda}$ is the vector of Lagrange multipliers, which
are widely used to enforce constraints in contact problems ([55, 56, 57, 58,
59, 60, 61, 62]).

System (35) contains the constraints on displacement level and is an index-3
system of differential-algebraic equations (DAE) [63]. The index of a DAE
system is the minimum number of differentiations needed to obtain ODEs for all
the unknowns. To obtain ODEs for the term ${\bf\Lambda}(t)={{\bf
L}^{\text{b}}(t)}^{{\mathsf{T}}}\boldsymbol{\lambda}$, equation (35b) has to
be differentiated once. A term $\dddot{\bf u}^{\text{b}}$ then appears in the
differentiated equation. Therefore, equation (35c) has to be differentiated
three times to solve for $\dddot{\bf u}^{\text{b}}$. If the constraint
equation (35c) is differentiated once or twice, the system can be reduced to
index-2 or index-1 DAE correspondingly. We next discuss numerical challenges
associated with solving an index-3 DAE system and with using straight finite
elements to model the bridge.

###  4.2 Numerical challenges

The solution of DAE systems usually requires careful consideration ([64, 65,
66, 67, 68]). When an index-3 DAE such as (35) is discretized using a time
integration scheme without numerical damping, spurious oscillations are
observed in wheel and bridge accelerations while using FEM to model the bridge
(Figures 6a-6d). This behavior, as well as oscillations in the Lagrange
multipliers, is due to the presence of the kinematic constraints. As explained
by Geradin et al.[69] and Cardona et al.[70], a weak instability in the second
time-derivative is observed when a standard Newmark scheme is used to
discretize a DAE system. There are two strategies to mitigate these
oscillations: 1) using numerical damping and 2) enforcing constraints at all
three levels (displacement, velocity, and acceleration), either explicitly or
through the correction steps (Section 6.3). The second approach is not
feasible in the case of the finite elements bridge model since standard finite
elements lack inter-element curvature continuity, and, as a result, the
acceleration constraint (51b) will be discontinuous. If the first approach is
employed, the introduction of numerical damping mitigates spurious
oscillations in the obtained response (Figures 6c-6d), similar to the effect
of the Bathe discretization scheme [71] in Reference [7]. While this approach
works well for mitigating the oscillations in the vertical direction, when a
curved bridge is introduced, the oscillations in the transverse direction are
amplified if straight FEM elements are used to approximate the curved shape of
the bridge model (Figure 7). Therefore, NURBS-based discretization of the
bridge is introduced in Section 5. Such discretization provides sufficient
inter-element continuity and also opens the path to enforcing acceleration
constraints explicitly (Section 6.2) and correcting wheel accelerations
through a simple constraint projection method (Section 6.3).

|
---|---
|
|

(a) Vertical bridge acceleration

(b) Vertical train accelerations

(c) Vertical bridge acceleration

(d) Vertical train accelerations

(e) Vertical bridge acceleration

(f) Vertical train accelerations

Figure 6:  Model: Simplified vehicle passing with speed $100\
$\mathrm{m}\mathrm{/}\mathrm{s}$$ over a single-span straight bridge modeled
with finite elements. Algorithmic feature tested: (a) and (b): Newmark scheme;
(c) and (d): The Generalized-$\alpha$ scheme without numerical damping
($\rho_{\infty}$= 1)444$\rho_{\infty}$ is the user-specified value of the
spectral radius in the high-frequency limit. The value $\rho_{\infty}=1$
corresponds to no dissipation, and the value $\rho_{\infty}=0$ results in
removing high-frequency response after the first time step (See Chung [72] and
Section 6.1).; (e) and (f): The Generalized-$\alpha$ scheme with slight
numerical damping ($\rho_{\infty}$= 0.9). Observations and Conclusions: The
presence of kinematic constraints causes unstable behavior when using
numerical schemes without numerical damping. The response can be stabilized by
introducing slight numerical damping.

|
---|---
|

(a) Vertical bridge acceleration (Eurocode EN 1990 (2002)[73] limits the value
of vertical bridge deck acceleration to $3.5\
$\mathrm{m}\mathrm{/}\mathrm{s}^{2}$$ for ballasted tracks)

(b) Vertical train accelerations

(c) Transverse contact force

(d) Transverse wheel acceleration

Figure 7:  Model: Simplified vehicle passing with speed $100\
$\mathrm{m}\mathrm{/}\mathrm{s}$$ over a 5-span bridge (Figure 3a) modeled
with finite elements. Algorithmic feature tested: The Generalized-$\alpha$
scheme (Section 6.1) with slight numerical damping ($\rho_{\infty}=0.9$).
Observations and Conclusions: Oscillations in the contact force (c) and wheel
accelerations (d) in the transverse direction are amplified due to the usage
of straight FEM elements in the bridge model.

##  5 NURBS-based beam formulation

In our past work on a two-dimensional VTSI algorithm [7], the Bathe method was
implemented to mitigate spurious oscillations caused by the lack of $C^{2}$
(curvature) continuity across beam elements in the bridge model. As explained
in Reference [52], the VTSI algorithm requires the second derivative of the
constraint matrix ${\bf L}^{\text{b}}$ to be twice continuously differentiable
spatially, for acceleration to be continuous in time. However, when the bridge
is modeled using standard finite elements, only $C^{1}$ continuity in the
nodes is provided.

In the present work, to provide sufficient inter-element continuity, an
isogeometric approach is employed for modeling the deck in the simplified
model and rails in the realistic model. This approach uses Non-Uniform
Rational B-Splines (NURBS) as a basis to describe the geometry and the
solution space [74]. Implementation of NURBS provides the desired $C^{2}$ (or
higher) continuity throughout the curve, as will be shown in Section 6. In
this study, the bridge is modeled using the three-dimensional isogeometric
beam formulation by Zhang et al. [75]. A thorough treatment of isogeometric
analysis in general can be found in Hughes et al.[74] and Cottrell et al.[76].
Algorithms for practical NURBS implementation are provided in Piegl et al.
[77]. For the use of NURBS in contact mechanics, see, for example, Agrawal et
al. [78]. Here, we provide a brief explanation of the modeling technique.

###  5.1 NURBS-based representation of geometry and solution fields

A NURBS curve can then be defined as (see Cottrell et al. [76] and Piegl et
al. [77])

${\bf\tilde{x}}(\xi)=\sum\limits_{i=1}^{n}R_{i,p}(\xi){\bf x}_{i}$ (36)

where ${\bf x}_{i}$ is the $i$-th control point, $n$ is the number of basis
functions, and $R_{i,p}(\xi)$ are the rational basis functions defined as

$R_{i,p}(\xi)=\frac{N_{i,p}(\xi)w_{i}}{\sum\limits_{j=1}^{n}N_{j,p}(\xi)w_{j}}=\frac{N_{i,p}(\xi)w_{i}}{W(\xi)}$
(37)

where $w_{i}$ are the control points’ weights and $N_{i,p}(\xi)$,
$i=1,2,...,n$, are the B-spline basis functions of degree $p$ defined on a
knot vector

$\Xi=\\{\underbrace{0,...,0}_{p+1},1,2,...,n^{e}-1,\underbrace{n^{e},...,n^{e}}_{p+1}\\}$
(38)

where the number of elements $n^{e}=n-p$ is equal to the number of nonzero
knot spans.

Equation 36 can be written in a matrix form as

${\bf\tilde{x}}(\xi)=\mathbf{R}^{\mathsf{T}}(\xi){\bf x}$ (39)

The geometry of the curve can either be represented exactly or, given a set of
data points, it can be approximated, for example, using least square
approximation. In either case, the NURBS basis is used to describe both
geometry and solution fields. Displacement and rotation fields can then be
approximated similarly to (39) as

$\begin{gathered}\tilde{\boldsymbol{\upsilon}}(\xi)=\sum\limits_{i=1}^{n}R_{i,p}(\xi)\boldsymbol{\upsilon}_{i}\\\
\tilde{\boldsymbol{\theta}}(\xi)=\sum\limits_{i=1}^{n}R_{i,p}(\xi)\boldsymbol{\theta}_{i}\end{gathered}$
(40)

where $\boldsymbol{\upsilon}$ and $\boldsymbol{\theta}$ are control
displacements and rotations.

Knots of the knot vector, $\Xi$, mesh the curve ${\bf\tilde{x}}$ into
elements. Each element ${\Omega}^{\text{e}}=[\xi_{i},\xi_{i+1}]$ is supported
by $p+1$ basis functions $R_{i-p,p},...,R_{i,p}$. The geometry, displacement,
and rotation fields at the point $\xi\in{\Omega}^{\text{e}}$ can then be
expressed as [75]

$\begin{gathered}{\bf\tilde{x}}(\xi)=\sum\limits_{k=i-p}^{i}R_{k,p}(\xi){\bf
x}_{k}={{{\bf R}^{\text{e}}}}^{\mathsf{T}}{\bf x}^{\text{e}}\\\
\tilde{\boldsymbol{\upsilon}}(\xi)=\sum\limits_{k=i-p}^{i}R_{k,p}(\xi)\boldsymbol{\upsilon}_{k}={{{\bf
R}^{\text{e}}}}^{\mathsf{T}}\boldsymbol{\upsilon}^{\text{e}}\\\
\tilde{\boldsymbol{\theta}}(\xi)=\sum\limits_{k=i-p}^{i}R_{k,p}(\xi)\boldsymbol{\theta}_{k}={{{\bf
R}^{\text{e}}}}^{\mathsf{T}}\boldsymbol{\theta}^{\text{e}}\\\ \end{gathered}$
(41)

###  5.2 Kinematics of the Timoshenko beam

Timoshenko beam theory is employed. Therefore, the generalized strains can be
expressed as [75, 79]

$\begin{gathered}\boldsymbol{\Theta}(s)=\left({\begin{array}[]{*{20}{c}}\boldsymbol{\varepsilon}(s)\\\
\boldsymbol{\beta}(s)\end{array}}\right)=\left({\begin{array}[]{c}\frac{d\boldsymbol{\upsilon}}{ds}-\boldsymbol{\theta}\times\textbf{t}\\\
\frac{d\boldsymbol{\theta}}{ds}\end{array}}\right)=\left({\begin{array}[]{l}\frac{d\upsilon_{t}}{ds}-\kappa\upsilon_{n}\\\\[2.84544pt]
\frac{d\upsilon_{n}}{ds}+\kappa\upsilon_{t}-\tau\upsilon_{b}-\theta_{b}\\\\[2.84544pt]
\frac{d\upsilon_{b}}{ds}+\tau\upsilon_{n}+\theta_{n}\\\\[2.84544pt]
\frac{d\theta_{t}}{ds}-\kappa\theta_{n}\\\\[2.84544pt]
\frac{d\theta_{n}}{ds}+\kappa\theta_{t}-\tau\theta_{b}\\\\[2.84544pt]
\frac{d\theta_{b}}{ds}-\tau\theta_{n}\end{array}}\right)\end{gathered}$ (42)

where $\boldsymbol{\varepsilon}$ is the vector of axial and the two shear
strains, $\boldsymbol{\beta}$ is the vector of one twisting and the two
bending strains, $\boldsymbol{\upsilon}$ and $\boldsymbol{\theta}$ are
independent displacement and rotation fields in Frenet basis
$\\{\textbf{t},\textbf{n},\textbf{b}\\}$, and $s$ is the arc-length parameter.
The parameter $\xi$ is mapped to the arc length as $\xi=\hat{\xi}(s)$, as
described in Section 2. Derivatives $\frac{d(\ )}{ds}$ are found through the
chain rule as $\frac{d(\ )}{ds}=\frac{d(\ )}{\xi}\frac{d\xi}{ds}=\frac{d(\
)}{\xi}\frac{1}{J}$, where the Jacobian is computed as
$J=\frac{ds}{d\xi}=\sqrt{(\frac{dx}{d\xi})^{2}+(\frac{dy}{d\xi})^{2}+(\frac{dz}{d\xi})^{2}}$.

###  5.3 Energy of the Timoshenko beam

The strain energy of the Timoshenko beam can be expressed as

$U=\frac{1}{2}\int\limits_{0}^{L}\left(EA\varepsilon_{t}^{2}+GA_{n}\varepsilon_{n}^{2}+GA_{b}\varepsilon_{b}^{2}+GI_{t}\beta_{t}^{2}+EI_{n}\beta_{n}^{2}+EI_{b}\beta_{b}^{2}\right)ds$
(43)

Having the generalized strain expressions (42) substituted into equation (43),
derivatives of (43) with respect to displacement and rotation fields are
found.

The kinetic energy of the beam is defined as

$T=\frac{1}{2}\int\limits_{0}^{L}\left(\begin{array}[]{l}\dot{\boldsymbol{\upsilon}}\\\
\dot{\boldsymbol{\theta}}\end{array}\right)^{T}\rho\left[\begin{array}[]{ll}A\textbf{I}&{}\\\
{}&\textbf{I}^{\text{b}}\end{array}\right]\left(\begin{array}[]{l}\dot{\boldsymbol{\upsilon}}\\\
\dot{\boldsymbol{\theta}}\end{array}\right)ds$ (44)

where $A$ is the cross-sectional area, I is an identity matrix, and
$\textbf{I}^{\text{b}}$ is a diagonal matrix with components
$\\{I_{t},I_{n},I_{b}\\}$, which are area moments of inertia about $t$, $n$,
and $b$ axes, correspondingly.

Work done by distributed loads, such as self-weight, can be found as

$V=\frac{1}{2}\int\limits_{0}^{L}\rho g\boldsymbol{\upsilon}_{b}ds$ (45)

The element load vector can then be found. Further details on the derivation
of member matrices can be found in Zhang et al. [75] and Fedorova [54].

##  6 Time integration

In this Section, three ways of solving the VTSI equations are described. The
three methods have their advantages and disadvantages. Some are easy to
implement in stand-alone codes. Others are better suited to be incorporated
into an existing structural analysis software.

###  6.1 Generalized-$\alpha$ method with numerical dissipation

As noted earlier, the system of VTSI equations (35) is an index-3 DAE system.
As the available research suggests, some amount of numerical damping is
necessary to avoid spurious high-frequency numerical oscillations in the
Lagrange multipliers in such systems [80]. Therefore, the Generalized-$\alpha$
method [72] is employed to discretize equations of motion in time. This scheme
provides accuracy at low-frequency and desired numerical damping at high-
frequency [80]. Other methods, such as Backward Difference Formulas and the
Bathe method [71] could also be used to this end; however, the
Generalized-$\alpha$ method is consistent with implementation in several
existing bridge analysis software, making it more convenient to integrate the
proposed VTSI algorithm, and is often used in discretization of constrained
systems ([81, 82]).

The VTSI equations of motion (35a)-(35b) are discretized in time using
averaging parameters $\alpha_{m}$ and $\alpha_{f}$ of the scheme as follows

$\displaystyle{\bf M}^{\text{t}}\ddot{\bf
u}^{\text{t}}_{{}_{n+1-\alpha_{m}}}+{\bf
C}^{\text{t}}({t}_{{}_{n+1-\alpha_{f}}})\dot{\bf
u}^{\text{t}}_{{}_{n+1-\alpha_{f}}}+{\bf
K}^{\text{t}}({t}_{{}_{n+1-\alpha_{f}}}){\bf
u}^{\text{t}}_{{}_{n+1-\alpha_{f}}}+{{\bf
L}^{\text{t}}}^{\mathsf{T}}\boldsymbol{\lambda}_{{}_{n+1-\alpha_{f}}}$
$\displaystyle={\bf P}^{\text{t}}({t}_{{}_{n+1-\alpha_{f}}})$ (46a)
$\displaystyle{\bf M}^{\text{b}}\ddot{\bf
u}^{\text{b}}_{{}_{n+1-\alpha_{m}}}+{\bf C}^{\text{b}}\dot{\bf
u}^{\text{b}}_{{}_{n+1-\alpha_{f}}}+{\bf K}^{\text{b}}{\bf
u}^{\text{b}}_{{}_{n+1-\alpha_{f}}}+\left({\bf
L}^{\text{b}}({t}_{{}_{n+1-\alpha_{f}}})\right)^{\mathsf{T}}\boldsymbol{\lambda}_{{}_{n+1-\alpha_{f}}}$
$\displaystyle={\bf P}^{\text{b}}$ (46b)

Even though Generalized-$\alpha$ is not an energy-preserving scheme, in order
to discretize the constraints, we assume that the work done by the constraint
forces over the time increment should vanish [69, 83] (as for an energy-
preserving scheme). This can be written using Equation (30) as

$W\approx(u_{i,\ n+1}-u_{i,\ n})\left(\frac{\partial C_{j}}{\partial
u_{i}}\right)_{{}_{n+1-\alpha_{f}}}\lambda_{j,\
{}_{n+1-\alpha_{f}}}\approx\left(C_{j,\ {}_{n+1}}-C_{j,\ n}\right)\lambda_{j,\
{}_{n+1-\alpha_{f}}}$ (47)

It follows from Equation (47) that enforcing constraint
${\left(C_{j}\right)}{}_{n+1}$ will ensure the vanishing of the work done by
constraint forces. Therefore, constraint equation (35c) is discretized in time
as

${\bf L}^{\text{t}}{\bf u}^{\text{t}}_{{}_{n+1}}+{\bf
L}^{\text{b}}(t_{n+1}){\bf u}^{\text{b}}_{{}_{n+1}}=\mathbf{0}$ (48)

The analysis steps are outlined in Procedure 1. The details of the time
discretization of Equations (46), as well as the notation of Procedure 1, can
be found in Fedorova [54]. Steps particular to 3D analysis and curved paths
are highlighted.

Figure 8 illustrates the results of VTSI analysis of the 5-span simplified
bridge (Figure 3a) modeled with NURBS and the simplified vehicle model (Figure
2). As can be seen in Figure 8c, implementation of NURBS helps to mitigate
spurious oscillations in the transverse direction (compare with Figure 7c,
where FEM elements are used). However, Figure 8b clearly demonstrates the
presence of acceleration discontinuity at time $t=0$. At this time, the
vehicle enters the bridge, and the curvature changes abruptly. The next
section discusses ways to avoid such discontinuity.

Procedure 1 Linear dynamic analysis of Vehicle-Track-Structure Interaction

1:Given: ${\bf P}^{\text{t}}$, ${\bf P}^{\text{b}}$, $\boldsymbol{\rho}(t)$,
$\Delta t$, $\alpha_{f}$, $\alpha_{m}$, $\beta$, $\gamma$, bridge geometry,
train geometry

2:Assemble ${\bf M}^{\text{t}}$, ${\bf L}^{\text{t}}$

3:Assemble ${\bf M}^{\text{b}}$, ${\bf C}^{\text{b}}$, ${\bf K}^{\text{b}}$

4:Assemble ${\bf{}\bar{M}}^{\text{b}}$

5:Compute ${\bf u}^{\text{t}}_{0}$ by linear static analysis. Set $\dot{\bf
u}^{\text{t}}_{0}=\bf{0}$ $\triangleright$ Train initial conditions

6:Compute ${\bf u}^{\text{b}}_{0}$ by linear static analysis. Set $\dot{\bf
u}^{\text{b}}_{0}=\bf{0}$ $\triangleright$ Bridge initial conditions

7:for $n\leftarrow 0$, NUMINC do

8: Compute ${t}_{n+1}=t_{n}+\Delta t$,
${t}_{{}_{n+1-\alpha_{f}}}=(1-{\alpha_{f}})t_{n+1}+{\alpha_{f}}t_{n}$

9: Compute the position of the wheel at the time step ${t}_{n+1}$,
${t}_{{}_{n+1-\alpha_{f}}}$

10: Compute Frenet frame properties ${\bf R}^{\text{F}}{}_{n+1-\alpha_{f}}$,
${\bf\ddot{x}}^{\text{F}}{}_{n+1-\alpha_{f}}$,
${\boldsymbol{\hat{\omega}}}^{\text{F}}{}_{n+1-\alpha_{f}}$,
${\boldsymbol{\dot{\hat{\omega}}}}^{\text{F}}{}_{n+1-\alpha_{f}}$

11: Obtain ${\bf L}^{\text{b}}({t}_{n+1})$, ${\bf
L}^{\text{b}}({t}_{{}_{n+1-\alpha_{f}}})$

12: Compute ${\bf a}^{\text{b}}{}_{n+1-\alpha_{f}}$

13: Compute ${\bf P}^{\text{t}}{}_{n+1-\alpha_{f}}$,
${\bf{}\bar{M}}^{\text{t}}{}_{n+1-\alpha_{f}}$, ${\bf
a}^{\text{t}}{}_{n+1-\alpha_{f}}$

14: Solve ${\bf{}\bar{M}}^{\text{t}}\left[{\bf\tilde{u}}^{\text{t}}\quad{\bf
A}^{\text{t}}\right]=\left[{\bf a}^{\text{t}}\quad{{\bf
L}^{\text{t}}}^{\mathsf{T}}\right]$

15: Solve ${\bf{}\bar{M}}^{\text{b}}\left[{\bf\tilde{u}}^{\text{b}}\quad{\bf
A}^{\text{b}}\right]=\left[{\bf a}^{\text{b}}\quad{\left({\bf
L}^{\text{b}}({t}_{{}_{n+1-\alpha_{f}}})\right)}^{{\mathsf{T}}}\right]$

16: Set ${\bf{}\bar{A}}=\left[{\bf L}^{\text{t}}{\bf A}^{\text{t}}+{\bf
L}^{\text{b}}({t}_{n+1}){\bf A}^{\text{b}}\right]$

17: Set $\boldsymbol{\bar{\rho}}={\bf
L}^{\text{t}}{\bf\tilde{u}}^{\text{t}}+{\bf
L}^{\text{b}}(t_{n+1}){\bf\tilde{u}}^{\text{b}}$

18: Solve
${\bf{}\bar{A}}\boldsymbol{\lambda}{}_{n+1-\alpha_{f}}=\boldsymbol{\bar{\rho}}$

19: Set ${\bf u}^{\text{t}}_{n+1}={\bf\tilde{u}}^{\text{t}}-{\bf
A}^{\text{t}}\boldsymbol{\lambda}{{}_{n+1-\alpha_{f}}}$; $\ddot{\bf
u}^{\text{t}}_{n+1}=\frac{1}{\Delta t^{2}\beta}\left({\bf
u}^{\text{t}}_{n+1}-{\bf u}^{\text{t}}_{n}-\Delta t\dot{\bf
u}^{\text{t}}_{n}\right)-\left(\frac{1}{2\beta}-1\right)\ddot{\bf
u}^{\text{t}}$;  $\dot{\bf u}^{\text{t}}_{n+1}=\dot{\bf
u}^{\text{t}}_{n}+\Delta t\left((1-\gamma)\ddot{\bf
u}^{\text{t}}_{n}+\gamma\ddot{\bf u}^{\text{t}}_{n+1}\right)$

20: Set ${\bf u}^{\text{b}}_{n+1}={\bf\tilde{u}}^{\text{b}}-{\bf
A}^{\text{b}}\boldsymbol{\lambda}{{}_{n+1-\alpha_{f}}}$; $\ddot{\bf
u}^{\text{b}}_{n+1}=\frac{1}{\Delta t^{2}\beta}\left({\bf
u}^{\text{b}}_{n+1}-{\bf u}^{\text{b}}_{n}-\Delta t\dot{\bf
u}^{\text{b}}_{n}\right)-\left(\frac{1}{2\beta}-1\right)\ddot{\bf
u}^{\text{b}}$;  $\dot{\bf u}^{\text{b}}_{n+1}=\dot{\bf
u}^{\text{b}}_{n}+\Delta t\left((1-\gamma)\ddot{\bf
u}^{\text{b}}_{n}+\gamma\ddot{\bf u}^{\text{b}}_{n+1}\right)$

21:end for

|
---|---
|

(a) Vertical bridge acceleration

(b) Vertical train accelerations

(c) Transverse contact force

(d) Transverse wheel acceleration

Figure 8:  Model: Simplified vehicle passing with speed $100\
$\mathrm{m}\mathrm{/}\mathrm{s}$$ over a 5-span bridge (Figure 3a) modeled
with NURBS of degree $p=3$. Algorithmic feature tested: Generalized-$\alpha$
scheme (Section 6.1) with slight numerical damping ($\rho_{\infty}=0.9$).
Observations and Conclusions: c) implementation of NURBS helps to mitigate
spurious oscillations in the transverse direction (compare with Figures 7c and
7d, where FEM elements are used. Note: The difference between the numerical
result and the theoretical is due to the fact that centripetal force is
calculated for the mass element only, not accounting for the effect of the
suspension element); b) acceleration discontinuity at time $t=0$: the vehicle
enters the bridge, and the curvature changes abruptly.

###  6.2 Explicit enforcement of acceleration constraints

One way to avoid the accumulation of error in obtained accelerations (see
Figures 6a-6b) is to enforce the acceleration constraint (51b) explicitly
instead of the displacement constraint (48). Such an approach is traditionally
used in multibody dynamics (for example, see Shabana [84]). In this case, the
VTSI equations of motion are the same as (35a)-(35b), and the acceleration
constraint can be discretized similarly to (48) as

${\bf L}^{\text{t}}\ddot{\bf
u}^{\text{t}}_{{}_{n+1}}+{\ddot{\bf{L}}}^{\text{b}}({t}_{n+1}){\bf
u}^{\text{b}}_{{}_{n+1}}+2\dot{\bf{L}}^{\text{b}}({t}_{n+1})\dot{\bf
u}^{\text{b}}_{{}_{n+1}}+{\bf L}^{\text{b}}({t}_{n+1})\ddot{\bf
u}^{\text{b}}_{{}_{n+1}}=0$ (49)

The system of Equations (46) and (49) can be rewritten similarly to how the
original system of Equations (46) and (48) was rewritten for the Procedure 1.

|
---|---
|

(a) Vertical bridge acceleration

(b) Vertical car accelerations

(c) Transverse contact force

(d) Transverse wheel accelerations

Figure 9:  Model: Simplified vehicle passing with speed $100\
$\mathrm{m}\mathrm{/}\mathrm{s}$$ over a 5-span bridge (Figure 3a) modeled
with NURBS of degree p = 3. Algorithmic feature tested: Acceleration
constraints are enforced explicitly (Equation (49)) at each time step. Wheel
velocity and acceleration are corrected only at $t=0$ using constraints (52).
Newmark scheme ($\alpha_{m}=0,\ {\alpha_{f}}=0$) is used. Observations and
Conclusions: Enforcement of acceleration constraint instead of displacement
constraint allows eliminating the spurious oscillations without employing the
Generalized-$\alpha$ scheme and by correcting only initial values of wheel
velocity and acceleration.

|
---|---
|

(a) Vertical bridge acceleration

(b) Vertical car accelerations

(c) Transverse contact force

(d) Transverse wheel accelerations

Figure 10:  Model: Simplified vehicle passing with speed $100\
$\mathrm{m}\mathrm{/}\mathrm{s}$$ over a 5-span bridge (Figure 3a) modeled
with NURBS degree p = 5. Algorithmic feature tested: Acceleration constraints
are enforced explicitly (Equation (49)) at each time step. Wheel velocity and
acceleration are corrected only at $t=0$ using constraints (52). Newmark
scheme ($\alpha_{m}=0,\ {\alpha_{f}}=0$) is used. Observations and
Conclusions: Enforcement of acceleration constraint instead of displacement
constraint allows eliminating the spurious oscillations without employing the
Generalized-$\alpha$ scheme and by correcting only initial values of wheel
velocity and acceleration. A higher degree of NURBS helps to further mitigate
the oscillations (compare with Figure 9).

As can be seen from Figures 9-10, enforcement of acceleration constraint
instead of displacement constraint allows to eliminate the spurious
oscillations without employing the Generalized-$\alpha$ scheme. Wheel velocity
and acceleration must be corrected only at $t=0$ using constraints (52).
Enforcing only acceleration constraints throughout the analysis leads to drift
in displacement constraints (see Figure 11a). Wheel displacements can be
corrected as

${\bf u}^{\text{w}}_{{}_{n+1}}={\bf L}^{\text{b}}(t_{n+1}){\bf
u}^{\text{b}}_{{}_{n+1}}$ (50)

Figure 11b shows the result of the correction at a single time step $t=0.75\
s$. In general, the correction of displacement constraints can be performed
periodically if needed.

|
---|---

(a) Displacement drift without correction

(b) Displacement drift with a correction of wheel displacement at $t=0.75\ s$

Figure 11:  Model: Simplified vehicle passing with speed $100\
$\mathrm{m}\mathrm{/}\mathrm{s}$$ over a 5-span bridge (Figure 3a) modeled
with NURBS degree $p=5$. Observations: Drift in the displacement constraint
occurs when only acceleration constraint (49) is enforced explicitly at each
time step.

###  6.3 Constraint projection

NURBS shape functions provide desired curvature continuity inside the bridge
span. However, at time $t=0$, when the train has not entered the bridge yet,
the continuity of velocity and acceleration constraints is not guaranteed.

Consider differentiating equation (35c) twice to obtain the constraints at
velocity and acceleration level

$\displaystyle{\bf L}^{\text{t}}\dot{\bf
u}^{\text{t}}+\dot{\bf{L}}^{\text{b}}{\bf u}^{\text{b}}+{\bf
L}^{\text{b}}\dot{\bf u}^{\text{b}}$ $\displaystyle=0$ (51a)
$\displaystyle{\bf L}^{\text{t}}\ddot{\bf
u}^{\text{t}}+{\ddot{\bf{L}}}^{\text{b}}{\bf
u}^{\text{b}}+2\dot{\bf{L}}^{\text{b}}\dot{\bf u}^{\text{b}}+{\bf
L}^{\text{b}}\ddot{\bf u}^{\text{b}}$ $\displaystyle=0$ (51b)

Recognizing that terms ${\bf L}^{\text{t}}\dot{\bf u}^{\text{t}}$ and ${\bf
L}^{\text{t}}\ddot{\bf u}^{\text{t}}$ correspond to wheel velocity and
acceleration, we can write

$\displaystyle\dot{\bf{u}}^{\text{w}}$
$\displaystyle=\dot{\bf{L}}^{\text{b}}{\bf u}^{\text{b}}+{\bf
L}^{\text{b}}\dot{\bf u}^{\text{b}}$ (52a)
$\displaystyle\ddot{\bf{u}}^{\text{w}}$
$\displaystyle={\ddot{\bf{L}}}^{\text{b}}{\bf
u}^{\text{b}}+2\dot{\bf{L}}^{\text{b}}\dot{\bf u}^{\text{b}}+{\bf
L}^{\text{b}}\ddot{\bf u}^{\text{b}}$ (52b)

At time $t=0$ wheel velocity and acceleration are equal to zero. At the same
time, the right-hand sides of equations (52) will not be equal to zero at this
time step, resulting in the discontinuities in wheel velocity and
acceleration. To overcome this, Equations (52) can be used to correct wheel
velocity and acceleration at time step $t=0$.

####  6.3.1 Simple example: Bridge modeled as a rigid cosine curve

Consider a simple example where the bridge in replaced with a rigid cosine
curve. By correcting the acceleration at time step $t=0$ using Equations (52),
we avoid discontinuities in the constraints at the acceleration level and
spurious oscillations in the accelerations (Figure 12).

|
---|---

(a) Vertical train accelerations without correction

(b) Vertical train accelerations with a correction (52b) of wheel acceleration
at $t=0$

Figure 12:  Model: Simulation of a vehicle running with speed $100\
$\mathrm{m}\mathrm{/}\mathrm{s}$$ over a rigid cosine curve. Algorithmic
feature tested: b) acceleration correction is performed only at $t=0$, while
displacement constraints are enforced explicitly (Equation (48)) at each time
step. Observations and Conclusions: By correcting the acceleration at time
step $t=0$ using Equations (52b), we avoid discontinuities in the constraints
at the acceleration level and spurious oscillations in the accelerations.

####  6.3.2 Example: Bridge modeled using NURBS

In a realistic case, when displacement constraints are enforced explicitly
using Equation (48), the correction only at time $t=0$ is not sufficient since
the error in accelerations accumulates and has to be corrected at each time
step using Equations (52). Figure 13 illustrates the case when velocity and
acceleration corrections are performed at each time step.

|
---|---
|

(a) Vertical bridge acceleration

(b) Vertical train accelerations

(c) Transverse contact force

(d) Transverse wheel acceleration

Figure 13:  Model: Simplified vehicle passing with speed $100\
$\mathrm{m}\mathrm{/}\mathrm{s}$$ over a 5-span bridge (Figure 3a) modeled
with NURBS of degree $p=3$. Algorithmic feature tested: Wheel velocity and
acceleration corrections (52) are performed at each time step, while
displacement constraints are enforced explicitly (Equation (48)). Newmark
scheme ($\alpha_{m}=0,\ {\alpha_{f}}=0$) is used. Observations and
Conclusions: The correction of velocity and acceleration constraints at each
time step helps to avoid the accumulation of errors in accelerations and makes
the use of the Newmark scheme possible.  Figure 14: Summary of algorithmic
choices

##  7 Numerical results

###  7.1 Simplified model

In this section, further numerical results for the simplified VTSI model
(Sections 3.1-3.2) are presented. These results are obtained using the
Generalized-$\alpha$ scheme with numerical dissipation (Section 6.1), since
among the presented choices, it is the most suitable way to integrate the VTSI
algorithm into an existing structural analysis software. Vertical bridge and
train accelerations, as well as transverse wheel acceleration and transverse
contact force, were already presented in Figure 8.

As can be seen from Figures 15 and 8, the VTSI analysis outputs such values as
displacement, accelerations, and contact forces (as well as velocities) for
all the model components. In practice, the obtained values have to be compared
with the code provisions.

|
---|---
|
|

(a) Vertical bridge displacements in the middle of the structure

(b) Transverse bridge displacements in the middle of the structure

(c) Vertical vehicle displacements

(d) Transverse vehicle displacements and rotations

(e) Vertical contact force

(f) Contact force in the direction $\theta_{n}$

Figure 15:  Model: Simplified vehicle passing with speed $100\
$\mathrm{m}\mathrm{/}\mathrm{s}$$ over a 5-span bridge (Figure 3a) modeled
with NURBS of degree $p=3$. Algorithmic feature tested: The
Generalized-$\alpha$ scheme (Section 6.1) with slight numerical damping
($\rho_{\infty}=0.9$).

###  7.2 Realistic model

A realistic train model from Section 3.3 and realistic bridge and track models
from Section 3.4 are used in this section. A one-car train (Figure 4) is
passing over a 15-span bridge (Figure 3b) with speed $100\
$\mathrm{m}\mathrm{/}\mathrm{s}$$, reaching the end of the bridge at time
$4.5\ $\mathrm{s}$$. The results in Figures 16 \- 22 are obtained using the
Generalized-$\alpha$ scheme with numerical dissipation. It was observed that
for models with many degrees of freedom, some numerical damping is needed,
even if other strategies, such as NURBS-based spatial discretization and
temporal constraint projection, are employed, since, in such models, the
oscillations arise not only due to the presence of the constraints but also
due to other high-frequency components in the solution. Therefore, if such a
method as constraint projection is used, a very small time step is required to
avoid the oscillations, which is not practical in realistic applications.

|
---|---

(a) Vertical bridge displacements

(b) Transverse bridge displacements

Figure 16:  Model: Realistic vehicle passing over a 15-span bridge (Figure 3b)
with speed $100\ $\mathrm{m}\mathrm{/}\mathrm{s}$$. Algorithmic feature
tested: The Generalized-$\alpha$ scheme (Section 6.1) with slight numerical
damping ($\rho_{\infty}=0.9$). Results: Bridge displacements at the middle of
the structure.

Eurocode EN 1990 (2002)[73] limits the value of vertical bridge deck
acceleration to $3.5\ $\mathrm{m}\mathrm{/}\mathrm{s}^{2}$$ (for a ballasted
track). The values presented in Figure 17a are well below the limit. The low
value can be attributed to the fact that a train composed of only a single car
was considered in this example to illustrate the modeling approaches described
in previous sections. A more significant response, including a resonance, is
expected when the bridge is analyzed under the trains with a standard number
of cars. To provide a ”very good” level of comfort to passengers, EN 1990[73]
limits vertical accelerations of a train car to $1\
$\mathrm{m}\mathrm{/}\mathrm{s}^{2}$$. The values in Figure 19a are below this
limit.

---

(a) Vertical bridge accelerations

(b) Transverse bridge accelerations

Figure 17:  Model: Realistic vehicle passing over a 15-span bridge (Figure 3b)
with speed $100\ $\mathrm{m}\mathrm{/}\mathrm{s}$$. Algorithmic feature
tested: The Generalized-$\alpha$ scheme (Section 6.1) with slight numerical
damping ($\rho_{\infty}=0.9$). Observations: The frequency of bridge
oscillations increases once realistic track and vehicle models are introduced
(compare, for example, with the simplified model’s results in Figure 13a).
However, the track provides additional damping to the system and redistributes
forces, resulting in smaller amplitude changes. The use of the
Generalized-$\alpha$ algorithm helps to mitigate spurious oscillations in this
realistic model with many degrees of freedom.

|
---|---

(a) Vertical car body displacements

(b) Transverse car body displacements

Figure 18:  Model: Realistic vehicle passing over a 15-span bridge (Figure 3b)
with speed $100\ $\mathrm{m}\mathrm{/}\mathrm{s}$$. Algorithmic feature
tested: The Generalized-$\alpha$ scheme (Section 6.1) with slight numerical
damping ($\rho_{\infty}=0.9$). Results: Train car body displacements.

|
---|---

(a) Vertical car body accelerations

(b) Transverse car body accelerations

Figure 19:  Model: Realistic vehicle passing over a 15-span bridge (Figure 3b)
with speed $100\ $\mathrm{m}\mathrm{/}\mathrm{s}$$. Algorithmic feature
tested: The Generalized-$\alpha$ scheme (Section 6.1) with slight numerical
damping ($\rho_{\infty}=0.9$). Results: Train car body accelerations.

---

(a) Vertical contact forces for the first wheelset

(b) Vertical accelerations of a wheelset contact point (right wheel of the
first wheelset)

Figure 20:  Model: Realistic vehicle passing over a 15-span bridge (Figure 3b)
with speed $100\ $\mathrm{m}\mathrm{/}\mathrm{s}$$. Algorithmic feature
tested: The Generalized-$\alpha$ scheme (Section 6.1) with slight numerical
damping ($\rho_{\infty}=0.9$). Observations: a) vertical contact force
increases at the right wheel and decreases at the left wheel of the first
wheelset while the train is passing over the curved spans; b) vertical
accelerations of a wheel contact point are obtained. Acceleration
discontinuity is observed when the vehicle enters the bridge, similar to the
simplified model (Figure 8b).

|
---|---

(a) Transverse contact forces for all the wheelset

(b) A sum of transverse contact forces and comparison with a theoretical value

(c) Transverse accelerations of a wheelset contact point (right wheel of the
first wheelset)

Figure 21:  Model: Realistic vehicle passing over a 15-span bridge (Figure 3b)
with speed $100\ $\mathrm{m}\mathrm{/}\mathrm{s}$$. Algorithmic feature
tested: The Generalized - $\alpha$ scheme (Section 6.1) with slight numerical
damping ($\rho_{\infty}=0.9$). Observations: Contact forces between wheels and
rails are obtained. a) transverse contact forces for all the wheelsets; b) sum
of contact forces resulting from the centrifugal acceleration of the train
compared with the theoretical value of the centripetal force (Note: The
difference between the numerical result and the theoretical is due to the fact
that centripetal force is calculated for a single mass element only, not
accounting for the effects of multiple bodies and suspensions); c) transverse
accelerations of a wheel contact point.

---

(a) Vertical contact force

(b) Transverse contact force

Figure 22:  Model: Realistic vehicle passing over a 15-span bridge (Figure 3b)
with speed $100\ $\mathrm{m}\mathrm{/}\mathrm{s}$$. Algorithmic feature
tested: The Generalized-$\alpha$ scheme (Section 6.1) with slight numerical
damping ($\rho_{\infty}=0.9$). Observations: Contact forces between the right
rail and the bridge are obtained in the middle of the structure. In the
present study, the connection between the bridge structure and the rails is
modeled using kinematic constraints. The use of The Generalized-$\alpha$
scheme and the rails formulated with NURBS allows for obtaining contact forces
without spurious oscillations.

##  8 Concluding remarks

An algorithm for three-dimensional dynamic vehicle-track-structure interaction
(VTSI) analysis has been discussed in this paper. It builds on our previous
work on a two-dimensional VTSI algorithm [7] and includes three main new
features: 1) a corotational approach to modeling vehicle dynamics, 2) NURBS-
based discretization of the vehicle path, 3) two alternative approaches to
solving an index-3 DAE system that are possible due to higher order continuity
provided by NURBS.

Employing the corotational approach, the train dynamics are conveniently
described with respect to Frenet frames along the train’s path on the bridge.
These Frenet frames serve as corotational frames. While displacements and
rotations of the train can be large, deformations within these corrotational
frames are assumed to be small, which allows the use of a linear (time-
varying) dynamic model. First, a simplified VTSI model is discussed in this
paper for clarity of demonstration, followed by an example of a realistic
model. Complete multibody train framework, as well as interacting bridge and
track models, are discussed in Reference [54].

While solving VTSI equations as an index-3 DAE system and employing standard
finite elements to model the track, spurious oscillations occur in obtained
contact forces and accelerations without appropriate algorithmic measures.
Moreover, if straight finite elements are used to approximate a curved bridge,
spurious oscillations are amplified in the transverse direction due to
curvature and rotation discontinuities. We employ isogeometric analysis to
model the bridge using NURBS. This helps mitigate oscillations in the
transverse direction. However, NURBS shape functions are twice continuously
differentiable only inside the bridge. This results in a discontinuity in
wheel accelerations at the beginning and end of the bridge. Correction steps
are proposed to avoid this discontinuity. Overall, three alternative
approaches to solving the VTSI system are discussed (see Figure 14 for the
summary).

The higher order continuity due to the use of NURBS allows for two alternative
approaches for solving the VTSI system. First, if the VTSI system is reduced
to index-1 DAE (when constraints are enforced at acceleration level), only
initial wheel velocity and acceleration need to be corrected to avoid the
discontinuity at the beginning of the bridge and the Newmark scheme can be
used for time integration (while avoiding numerical damping). This correction
ensures consistent initial conditions at all three levels (displacement,
velocity, and acceleration) of constraints. When displacement constraints are
not enforced explicitly, drift in displacement compatibility develops. The
observed drift is not significant, but a periodic wheel displacement
correction based on the kinematic constraint can be used as needed. Second, it
has also been demonstrated that index-3 VTSI equations can be solved using the
Newmark scheme if wheel velocities and accelerations are corrected at each
time step using a simple non-iterative constraint projection method.

The proposed approach allows for a convenient description of vehicle dynamics
through corotational frames and for solving three-dimensional VTSI as an
index-3 DAE system with or without numerical damping. The algorithm can be
applied to other multibody systems (with large displacements and rotations,
but small deformations) coupled to deformable structures by time-varying
kinematic constraints. The use of kinematic constraints leads to the
modularity of the proposed VTSI algorithm and its straightforward integration
into existing finite element software.

## Acknowledgments

The first author gratefully acknowledges a Ph.D. scholarship from LARSA Inc.,
Melville, NY.

##  Appendix A Simplified vehicle matrices

This appendix includes matrices for a simplified vehicle model presented in
Section 3.1. Here, ${\bf M}^{\text{t}}$, ${\bf C}^{\text{t}}(t)$, and ${\bf
K}^{\text{t}}(t)$ are the mass, damping, and stiffness matrices. ${\bf
L}^{\text{t}}$ is the kinematic constraint matrix used to constrain wheels to
the rails. $\boldsymbol{\lambda}$ is a vector of constraint forces in
transverse, vertical, and rolling directions. ${\bf P}^{\text{t}}$ is the load
vector.

${\bf
M}^{\text{t}}=\left[{\begin{array}[]{*{20}{c}}{m^{\text{w}}+m^{\text{c}}}&{}\hfil&{{-m^{\text{c}}}{l_{0}}}&{}\hfil\\\
{}\hfil&{{m^{\text{w}}}}&{}\hfil&{}\hfil\\\
{{-m^{\text{c}}}{l_{0}}}&{}\hfil&{{m^{\text{c}}}{l_{0}}^{2}+{I^{\text{w}}}+{I^{\text{c}}}}&{}\hfil\\\
{}\hfil&{}\hfil&{}\hfil&{{m^{\text{c}}}}\end{array}}\right]$ (A.1)

${\bf C}^{\text{t}}=2m^{\text{w}}{{\bf
T}^{\text{w}}}^{\mathsf{T}}{\boldsymbol{\hat{\omega}}}^{\text{F}}{\bf
T}^{\text{w}}+2m^{\text{c}}{{\bf
T}^{\text{c}}}^{\mathsf{T}}{\boldsymbol{\hat{\omega}}}^{\text{F}}{\bf
T}^{\text{c}}=\left[{\begin{array}[]{*{20}{c}}{}\hfil&{-2{m^{\text{w}}}{\omega}^{\text{F}}_{1}}&{}\hfil&{-2{m^{\text{c}}}{\omega}^{\text{F}}_{1}}\\\
{2{m^{\text{w}}}{\omega}^{\text{F}}_{1}}&{}\hfil&{}\hfil&{}\hfil\\\
{}\hfil&{}\hfil&{}\hfil&{2{m^{\text{c}}}{\omega}^{\text{F}}_{1}{l_{0}}}\\\
{2{m^{\text{c}}}{\omega}^{\text{F}}_{1}}&{}\hfil&{-2{m^{\text{c}}}{\omega}^{\text{F}}_{1}{l_{0}}}&{}\hfil\end{array}}\right]$
(A.2)

$\begin{gathered}{\bf K}^{\text{t}}=m^{\text{w}}{{\bf
T}^{\text{w}}}^{\mathsf{T}}{\boldsymbol{\dot{\hat{\omega}}}}^{\text{F}}{\bf
T}^{\text{w}}+m^{\text{c}}{{\bf
T}^{\text{c}}}^{\mathsf{T}}{\boldsymbol{\dot{\hat{\omega}}}}^{\text{F}}{\bf
T}^{\text{c}}+m^{\text{w}}{{\bf
T}^{\text{w}}}^{\mathsf{T}}{{\boldsymbol{\hat{\omega}}}^{\text{F}}}^{2}{\bf
T}^{\text{w}}+m^{\text{c}}{{\bf
T}^{\text{c}}}^{\mathsf{T}}{{\boldsymbol{\hat{\omega}}}^{\text{F}}}^{2}{\bf
T}^{\text{c}}=\\\
\left[{\begin{array}[]{*{20}{c}}{-\left({{m^{\text{w}}}+{m^{\text{c}}}}\right)\left({{\omega}^{\text{F}}_{1}}^{2}+{{\omega}^{\text{F}}_{3}}^{2}\right)}&{-{m^{\text{w}}}\left({\dot{\omega}}^{\text{F}}_{1}-{\omega}^{\text{F}}_{2}{\omega}^{\text{F}}_{3}\right)}&{{m^{\text{c}}}\left({{\omega}^{\text{F}}_{1}}^{2}+{{\omega}^{\text{F}}_{3}}^{2}\right){l_{0}}}&{-{m^{\text{c}}}\left({\dot{\omega}}^{\text{F}}_{1}-{\omega}^{\text{F}}_{2}{\omega}^{\text{F}}_{3}\right)}\\\
{{m^{\text{w}}}\left({{\dot{\omega}}^{\text{F}}_{1}+{\omega}^{\text{F}}_{2}{\omega}^{\text{F}}_{3}}\right)}&{-{m^{\text{w}}}\left({{\omega}^{\text{F}}_{1}}^{2}+{{\omega}^{\text{F}}_{2}}^{2}\right)+{k^{\text{s}}}}&0&{-{k^{\text{s}}}}\\\
{{m^{\text{c}}}\left({{\omega}^{\text{F}}_{1}}^{2}+{{\omega}^{\text{F}}_{3}}^{2}\right){l_{0}}}&0&{-{m^{\text{c}}}\left({{\omega}^{\text{F}}_{1}}^{2}+{{\omega}^{\text{F}}_{3}}^{2}\right){l_{0}}^{2}}&{{m^{\text{c}}}\left({{\dot{\omega}}^{\text{F}}_{1}-{\omega}^{\text{F}}_{2}{\omega}^{\text{F}}_{3}}\right){l_{0}}}\\\
{{m^{\text{c}}}\left({{\dot{\omega}}^{\text{F}}_{1}+{\omega}^{\text{F}}_{2}{\omega}^{\text{F}}_{3}}\right)}&{-{k^{\text{s}}}}&{-{m^{\text{c}}}\left({{\dot{\omega}}^{\text{F}}_{1}+{\omega}^{\text{F}}_{2}{\omega}^{\text{F}}_{3}}\right){l_{0}}}&{-{m^{\text{c}}}\left({{\omega}^{\text{F}}_{1}}^{2}+{{\omega}^{\text{F}}_{2}}^{2}\right)+{k^{\text{s}}}}\end{array}}\right]\end{gathered}$
(A.3)

${\bf L}^{\text{t}}=\left[{\begin{array}[]{*{20}{c}}{-1}&0&0\\\ 0&{-1}&0\\\
0&0&{-1}\\\ 0&0&0\end{array}}\right]$ (A.4)

$\boldsymbol{\lambda}=\left({\begin{array}[]{*{20}{c}}{{\lambda_{y}}}\\\
{{\lambda_{z}}}\\\ {{\lambda_{{\theta_{x}}}}}\end{array}}\right)$ (A.5)

$\begin{gathered}{\bf P}^{\text{t}}=-\left(m^{\text{w}}{{\bf
T}^{\text{w}}}^{\mathsf{T}}+m^{\text{c}}{{\bf
T}^{\text{c}}}^{\mathsf{T}}\right){{\bf
R}^{\text{F}}}^{\mathsf{T}}{\bf\ddot{x}}^{\text{F}}-m^{\text{c}}{{\bf
T}^{\text{c}}}^{\mathsf{T}}{\boldsymbol{\dot{\hat{\omega}}}}^{\text{F}}\left({\begin{array}[]{*{20}{c}}0\\\
0\\\ l_{0}\end{array}}\right)-m^{\text{c}}{{\bf
T}^{\text{c}}}^{\mathsf{T}}{{\boldsymbol{\hat{\omega}}}^{\text{F}}}^{2}\left({\begin{array}[]{*{20}{c}}0\\\
0\\\ l_{0}\end{array}}\right)\\\ -(m^{\text{w}}g{{\bf
T}^{\text{w}}}^{\mathsf{T}}+m^{\text{c}}g{{\bf
T}^{\text{c}}}^{\mathsf{T}})({{\bf R}^{\text{F}}(t)}^{\mathsf{T}}-{{\bf
R}^{\text{F}}(0)}^{\mathsf{T}})\left({\begin{array}[]{*{20}{c}}0\\\ 0\\\
1\end{array}}\right)\end{gathered}$ (A.6)

The vehicle matrices listed above are used in the simplified vehicle equations
of motion (32). Complete derivations of the equations of motion can be found
in Fedorova [52].

## References

  * Bowers [2013] M. Bowers, California High-Speed Train Project. Technical Memorandum. Track-Structure Interaction TM 2.10.10, Technical Report, Parsons Brinckerhoff, California High-Speed Rail Authority, 2013.
  * Eur [2003] EN 1991-2. Eurocode 1: Actions on structures - Part 2: Traffic Loads on bridges, Technical Report, European Committee for Standardization, 2003.
  * Fedorova et al. [2019] M. Fedorova, M. V. Sivaselvan, J. Tauberer,  On the right track,  Bridge design & engineering 94 (2019) 72–73.
  * Yang et al. [2004] Y.-B. Yang, J. Yau, Y. Wu, Vehicle-Bridge Interaction Dynamics With Applications to High-Speed Railways, World Scientific, 2004.
  * Shabana et al. [2007] A. A. Shabana, K. E. Zaazaa, H. Sugiyama, Railroad Vehicle Dynamics: A Computational Approach, CRC Press, 2007. URL: https://www.taylorfrancis.com/books/9781420045857.
  * Zhai et al. [2019] W. Zhai, Z. Han, Z. Chen, L. Ling, S. Zhu,  Train-track-bridge dynamic interaction: a state-of-the-art review,  Vehicle System Dynamics 57 (2019) 984–1027. doi:10.1080/00423114.2019.1605085.
  * Fedorova and Sivaselvan [2017] M. Fedorova, M. V. Sivaselvan,  An algorithm for dynamic vehicle-track-structure interaction analysis for high-speed trains,  Engineering Structures 148 (2017) 857–877. doi:10.1016/j.engstruct.2017.05.065.
  * Neves et al. [2012] S. Neves, A. Azevedo, R. Calcada,  A direct method for analyzing the vertical vehicle-structure interaction,  Engineering Structures 34 (2012) 414–420.
  * Xia and Zhang [2005] H. Xia, N. Zhang,  Dynamic analysis of railway bridge under high-speed trains,  Computers & Structures 83 (2005) 1891–1901.
  * Dimitrakopoulos and Zeng [2015] E. Dimitrakopoulos, Q. Zeng,  A three-dimensional dynamic analysis scheme for the interaction between trains and curved railway bridges,  Computers and Structures 149 (2015) 43–60.
  * Salcher and Adam [2015] P. Salcher, C. Adam,  Modeling of dynamic train-bridge interaction in high-speed railways,  Acta Mechanica 226 (2015) 2473–2495.
  * Arvidsson et al. [2019] T. Arvidsson, A. Andersson, R. Karoumi,  Train running safety on non-ballasted bridges,  International Journal of Rail Transportation 7 (2019) 1–22. doi:10.1080/23248378.2018.1503975.
  * Bettinelli et al. [2023] L. Bettinelli, A. Stollwitzer, J. Fink,  Numerical Study on the Influence of Coupling Beam Modeling on Structural Accelerations during High-Speed Train Crossings,  Applied Sciences (Switzerland) 13 (2023). doi:10.3390/app13158746.
  * Glatz and Fink [2021] B. Glatz, J. Fink,  A redesigned approach to the additional damping method in the dynamic analysis of simply supported railway bridges,  Engineering Structures 241 (2021). doi:10.1016/j.engstruct.2021.112415.
  * Greco and Lonetti [2018] F. Greco, P. Lonetti,  Numerical formulation based on moving mesh method for vehicle-bridge interaction,  Advances in Engineering Software 121 (2018) 75–83. doi:10.1016/j.advengsoft.2018.03.013.
  * Mosayebi et al. [2021] S. A. Mosayebi, J. A. Zakeri, M. Esmaeili,  Investigation and comparison of dynamic interaction models of vehicle-track systems,  Journal of Theoretical and Applied Mechanics (Poland) 58 (2021) 273–284. doi:10.15632/jtam-pl/115971.
  * Gou et al. [2018] H. Gou, X. Shi, W. Zhou, K. Cui, Q. Pu,  Dynamic performance of continuous railway bridges: Numerical analyses and field tests,  Proceedings of the Institution of Mechanical Engineers, Part F: Journal of Rail and Rapid Transit 232 (2018) 936–955. doi:10.1177/0954409717702019.
  * Xiao et al. [2020] Y. Xiao, X. Luo, J. Liu, K. Wang,  Dynamic Response of Railway Bridges under Heavy-Haul Freight Trains,  Advances in Civil Engineering 2020 (2020). doi:10.1155/2020/7486904.
  * Xu and Zhai [2019] L. Xu, W. Zhai,  A three-dimensional model for train-track-bridge dynamic interactions with hypothesis of wheel-rail rigid contact,  Mechanical Systems and Signal Processing 132 (2019) 471–489. doi:10.1016/j.ymssp.2019.04.025.
  * Xiao and Ren [2019] X. Xiao, W. X. Ren,  A Versatile 3D Vehicle-Track-Bridge Element for Dynamic Analysis of the Railway Bridges under Moving Train Loads,  International Journal of Structural Stability and Dynamics 19 (2019). doi:10.1142/S0219455419500500.
  * Eroğlu et al. [2022] M. Eroğlu, M. A. Koç, İ. Esen, R. Kozan,  Train-structure interaction for high-speed trains using a full 3D train model,  Journal of the Brazilian Society of Mechanical Sciences and Engineering 44 (2022). doi:10.1007/s40430-021-03338-1.
  * Zhou and Zhang [2023] Z.-j. Zhou, N. Zhang,  Statistical study on the derailment factor considering the weak correlation of lateral-vertical wheel-rail forces of a vehicle-bridge system,  Vehicle System Dynamics (2023). doi:10.1080/00423114.2023.2170810.
  * Yu and Mao [2018] Z.-w. Yu, J.-f. Mao,  A stochastic dynamic model of train-track-bridge coupled system based on probability density evolution method,  Applied Mathematical Modelling 59 (2018) 205–232. doi:10.1016/j.apm.2018.01.038.
  * He et al. [2020] Q. He, C. Cai, S. Zhu, M. Wang, K. Wang, W. Zhai,  Key parameter selection of suspended monorail system based on vehicle-bridge dynamical interaction analysis,  Vehicle System Dynamics 58 (2020) 339–356. doi:10.1080/00423114.2019.1577470.
  * Wang et al. [2019] Z. L. Wang, Y. L. Xu, G. Q. Li, Y. B. Yang, S. W. Chen, X. L. Zhang,  Modelling and validation of coupled high-speed maglev train-and-viaduct systems considering support flexibility,  Vehicle System Dynamics 57 (2019) 161–191. doi:10.1080/00423114.2018.1450517.
  * Li et al. [2022] M. Li, S. Luo, W. Ma, T. Li, D. Gao, Z. Xu,  Experimental and numerical investigations of the dynamic responses of low and medium speed maglev train-track-bridge coupled system,  Vehicle System Dynamics 60 (2022) 1555–1578. doi:10.1080/00423114.2020.1864417.
  * Zhang et al. [2023] M. Zhang, C. Yuan, W. Ma, S. Luo,  Effect of suspension form on the vehicle-bridge coupled vibration of the maglev vehicle,  Vehicle System Dynamics (2023). doi:10.1080/00423114.2023.2180756.
  * Shi et al. [2022] Y. Shi, W. Ma, M. Li, Z. Xu,  Research on dynamics of a new high-speed maglev vehicle,  Vehicle System Dynamics 60 (2022) 721–742. doi:10.1080/00423114.2020.1838568.
  * Chen et al. [2021] Z. Chen, L. Bi, J. Zhao,  Comparison of single-pier settlement model and multi-pier settlement model in solving train-track-bridge interaction,  Vehicle System Dynamics 59 (2021) 1484–1508. doi:10.1080/00423114.2020.1763406.
  * Chen et al. [2022] L. K. Chen, A. Kurtulus, Y. F. Dong, E. Taciroglu, L. Z. Jiang,  Velocity pulse effects of near-fault earthquakes on a high-speed railway vehicle-ballastless track-benchmark bridge system,  Vehicle System Dynamics 60 (2022) 2963–2987. doi:10.1080/00423114.2021.1933546.
  * Montenegro et al. [2020] P. A. Montenegro, R. Heleno, H. Carvalho, R. Calçada, C. J. Baker,  A comparative study on the running safety of trains subjected to crosswinds simulated with different wind models,  Journal of Wind Engineering and Industrial Aerodynamics 207 (2020). doi:10.1016/j.jweia.2020.104398.
  * Xin et al. [2020] L. Xin, X. Li, Y. Zhu, M. Liu,  Uncertainty and sensitivity analysis for train-ballasted track-bridge system,  Vehicle System Dynamics 58 (2020) 453–471. doi:10.1080/00423114.2019.1584678.
  * Su et al. [2022] C. Su, Z. Wu, J. Xian,  Stochastic dynamic analysis of vehicle-bridge coupled systems with nonlinear Hertz contacts by explicit time-domain method,  Vehicle System Dynamics 60 (2022) 1579–1601. doi:10.1080/00423114.2020.1864418.
  * Diana and Cheli [1989] G. Diana, F. Cheli,  Dynamic Interaction of Railway Systems with Large Bridges,  Vehicle System Dynamics 18 (1989) 71–106.
  * Kwark et al. [2004] J. W. Kwark, E. S. Choi, Y. J. Kim, B. S. Kim, S. I. Kim,  Dynamic behavior of two-span continuous concrete bridges under moving high-speed train,  Computers & Structures 82 (2004) 463–474.
  * Zhang and Xia [2013] N. Zhang, H. Xia,  Dynamic analysis of coupled vehicle-bridge system based on inter-system iteration method,  Computers & Structures 114-115 (2013) 26–34.
  * Ticona Melo et al. [2018] L. R. Ticona Melo, T. N. Bittencourt, D. Ribeiro, R. Calçada,  Dynamic Response of a Railway Bridge to Heavy Axle-Load Trains Considering Vehicle-Bridge Interaction,  International Journal of Structural Stability and Dynamics 18 (2018). doi:10.1142/S0219455418500104.
  * Ticona Melo et al. [2020] L. R. Ticona Melo, D. Ribeiro, R. Calçada, T. N. Bittencourt,  Validation of a vertical train-track-bridge dynamic interaction model based on limited experimental data,  Structure and Infrastructure Engineering 16 (2020) 181–201. doi:10.1080/15732479.2019.1605394.
  * Xu et al. [2022] L. Xu, Z. Li, Y. Zhao, Z. Yu, K. Wang,  Modelling of vehicle-track related dynamics: a development of multi-finite-element coupling method and multi-time-step solution method,  Vehicle System Dynamics 60 (2022) 1097–1124. doi:10.1080/00423114.2020.1847298.
  * Lu et al. [2020] X. Lu, C. W. Kim, K. C. Chang,  Finite Element Analysis Framework for Dynamic Vehicle-Bridge Interaction System Based on ABAQUS,  International Journal of Structural Stability and Dynamics 20 (2020). doi:10.1142/S0219455420500340.
  * Liu et al. [2022] L. Liu, Z. Zuo, J. D. Yau, S. Urushadze,  A simplified method to assess vehicle-bridge interaction for train-induced vibration of light-weight railway bridges,  Journal of the Chinese Institute of Engineers, Transactions of the Chinese Institute of Engineers,Series A 45 (2022) 651–660. doi:10.1080/02533839.2022.2126400.
  * Li et al. [2021] H. Li, T. Wang, G. Wu,  Dynamic response prediction of vehicle-bridge interaction system using feedforward neural network and deep long short-term memory network,  Structures 34 (2021) 2415–2431. doi:10.1016/j.istruc.2021.09.008.
  * Zhou et al. [2023] J. Zhou, C. Huang, J. Deng, J. Zhang, L. Zhang,  A co-simulation approach for straddle monorail vehicle-bridge interaction subjected to nonlinear excitation,  Advances in Engineering Software 180 (2023). doi:10.1016/j.advengsoft.2023.103458.
  * Stefanidou and Paraskevopoulos [2022] S. P. Stefanidou, E. A. Paraskevopoulos,  Seismic fragility analysis of railway reinforced concrete bridges considering real-time vehicle-bridge interaction with the aid of co-simulation techniques,  Earthquake Engineering and Structural Dynamics 51 (2022) 2137–2161. doi:10.1002/eqe.3657.
  * Sivaselvan et al. [2014] M. V. Sivaselvan, J. Tauberer, A. Karakaplan,  Dynamic vehicle-track-structure interaction analysis using lagrange multipliers,  Proceedings of the Istanbul Bridge Conference, 2014.
  * Zhu et al. [2015] D. Y. Zhu, Y. H. Zhang, H. Ouyang,  A linear complementarity method for dynamic analysis of bridges under moving vehicles considering separation and surface roughness,  Computers & Structures 154 (2015) 135–144.
  * Felippa and Haugen [2005] C. A. Felippa, B. Haugen,  A unified formulation of small-strain corotational finite elements: I. Theory,  Computer Methods in Applied Mechanics and Engineering 194 (2005) 2285–2335. doi:10.1016/j.cma.2004.07.035.
  * You et al. [2024] P. You, Z. Liu, Z. Ma,  A contact formulation using local frame of SE(3) group for corotational beam against rigid solid via moving signed distance field,  Computer Methods in Applied Mechanics and Engineering 418 (2024). doi:10.1016/j.cma.2023.116532.
  * Farin [1997] G. E. Farin, Curves and surfaces for computer-aided geometric design : a practical guide, Academic Press, 1997.
  * Tapp [2016] K. Tapp, Differential Geometry of Curves and Surfaces, Springer International Publishing, 2016.
  * Murray et al. [1994] R. Murray, Z. Li, S. Sastry, A mathematical introduction to robotic manipulation, CRC press, 1994.
  * Fedorova [2017] M. Fedorova, An algorithm for dynamic vehicle-track-structure interaction analysis for high-speed trains, Master’s thesis, State University of New York at Buffalo, 2017. https://search.proquest.com/docview/1877654821.
  * Harris [2009] G. Harris, California High-Speed Train Project. Alignment Design Standards for High-Speed Train Operation. TM 2.1.2, Technical Report, California High-Speed Rail Authority, 2009.
  * Fedorova [2024] M. Fedorova, A modular algorithm for three-dimensional dynamic vehicle-track-structure interaction analysis for high-speed trains using Lagrange multipliers, Ph.D. Dissertation, State University of New York at Buffalo, Buffalo, NY, 2024.
  * Kothari and Krause [2022] H. Kothari, R. Krause,  A generalized multigrid method for solving contact problems in Lagrange multiplier based unfitted Finite Element method,  Computer Methods in Applied Mechanics and Engineering 392 (2022). doi:10.1016/j.cma.2022.114630. arXiv:2106.16107.
  * Duong et al. [2017] T. X. Duong, F. Roohbakhshan, R. A. Sauer,  A new rotation-free isogeometric thin shell formulation and a corresponding continuity constraint for patch boundaries,  Computer Methods in Applied Mechanics and Engineering (2017). doi:10.1016/j.cma.2016.04.008.
  * Meier et al. [2016] C. Meier, A. Popp, W. A. Wall,  A finite element approach for the line-to-line contact interaction of thin beams with arbitrary orientation,  Computer Methods in Applied Mechanics and Engineering 308 (2016) 377–413. doi:10.1016/j.cma.2016.05.012. arXiv:1603.09227.
  * Kadapa et al. [2016] C. Kadapa, W. G. Dettmer, D. Perić,  A fictitious domain/distributed Lagrange multiplier based fluid-structure interaction scheme with hierarchical B-Spline grids,  Computer Methods in Applied Mechanics and Engineering 301 (2016) 1–27. doi:10.1016/j.cma.2015.12.023.
  * Matzen and Bischoff [2016] M. E. Matzen, M. Bischoff,  A weighted point-based formulation for isogeometric contact,  Computer Methods in Applied Mechanics and Engineering 308 (2016) 73–95. doi:10.1016/j.cma.2016.04.010.
  * Zhang et al. [2011] H. W. Zhang, J. Y. Li, S. H. Pan,  New second-order cone linear complementarity formulation and semi-smooth Newton algorithm for finite element analysis of 3D frictional contact problem,  Computer Methods in Applied Mechanics and Engineering 200 (2011) 77–88. doi:10.1016/j.cma.2010.07.014.
  * Simeon [2006] B. Simeon,  On Lagrange multipliers in flexible multibody dynamics,  Computer Methods in Applied Mechanics and Engineering 195 (2006) 6993–7005. doi:10.1016/j.cma.2005.04.015.
  * Baumgarte [1972] J. Baumgarte,  Stabilization of constraints and integrals of motion in dynamical systems,  Computer Methods in Applied Mechanics and Engineering 1 (1972) 1–16. doi:10.1016/0045-7825(72)90018-7.
  * Brenan et al. [1996] K. E. Brenan, S. L. V. Campbell, S. L. Campbell, L. R. Petzold, Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations, Classics in Applied Mathematics, Society for Industrial and Applied Mathematics, 1996.
  * Matthies et al. [2006] H. G. Matthies, R. Niekamp, J. Steindorf,  Algorithms for strong coupling procedures,  Computer Methods in Applied Mechanics and Engineering 195 (2006) 2028–2049. doi:10.1016/j.cma.2004.11.032.
  * Arnold et al. [2006] M. Arnold, A. Fuchs, C. Führer,  Efficient corrector iteration for DAE time integration in multibody dynamics,  Computer Methods in Applied Mechanics and Engineering 195 (2006) 6958–6973. doi:10.1016/j.cma.2005.02.034.
  * Borri et al. [2006] M. Borri, L. Trainelli, A. Croce,  The Embedded Projection Method: A general index reduction procedure for constrained system dynamics,  Computer Methods in Applied Mechanics and Engineering 195 (2006) 6974–6992. doi:10.1016/j.cma.2005.03.010.
  * Betsch [2005] P. Betsch,  The discrete null space method for the energy consistent integration of constrained mechanical systems. Part I: Holonomic constraints,  Computer Methods in Applied Mechanics and Engineering (2005). doi:10.1016/j.cma.2005.01.004. arXiv:arXiv:1011.1669v3.
  * Fisette and Vaneghem [1996] P. Fisette, B. Vaneghem,  Numerical integration of multibody system dynamic equations using the coordinate partitioning method in an implicit Newmark scheme,  Computer Methods in Applied Mechanics and Engineering (1996). doi:10.1016/0045-7825(95)00926-4.
  * Géradin and Cardona [2001] M. Géradin, A. Cardona, Flexible multibody dynamics: a finite element approach, John Wiley, 2001.
  * Cardona and Geradin [1989] A. Cardona, M. Geradin,  Time integration of the equations of motion in mechanism analysis,  Computers & Structures 33 (1989) 801–820. doi:10.1016/0045-7949(89)90255-1.
  * Bathe and Baig [2005] K. Bathe, M. Baig,  On a composite implicit time integration procedure for nonlinear dynamics,  Computers & Structures 83 (2005) 2513–2524.
  * Chung and Hulbert [1994] J. Chung, G. M. Hulbert,  A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-$\alpha$ method,  Journal of Applied Mechanics 60 (1994) 371–375.
  * CEN/TC 250 [2005] CEN/TC 250, EN 1990 (2002): Eurocode - Basis of structural design, 2005.
  * Hughes et al. [2005] T. J. Hughes, J. A. Cottrell, Y. Bazilevs,  Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement,  Computer Methods in Applied Mechanics and Engineering 194 (2005) 4135–4195. doi:10.1016/j.cma.2004.10.008.
  * Zhang et al. [2016] G. Zhang, R. Alberdi, K. Khandelwal,  Analysis of three-dimensional curved beams using isogeometric approach,  Engineering Structures 117 (2016) 560–574.
  * Cottrell et al. [2009] J. A. Cottrell, T. J. Hughes, Y. Bazilevs, Isogeometric Analysis: Toward Integration of CAD and FEA, Wiley, 2009.
  * Piegl and Tiller [1997] L. Piegl, W. Tiller, The NURBS Book (2Nd Ed.), Springer-Verlag New York, Inc., New York, NY, USA, 1997.
  * Agrawal and Gautam [2020] V. Agrawal, S. S. Gautam,  Varying-order NURBS discretization: An accurate and efficient method for isogeometric analysis of large deformation contact problems,  Computer Methods in Applied Mechanics and Engineering 367 (2020). doi:10.1016/j.cma.2020.113125. arXiv:1903.05859.
  * Tabarrok et al. [1988] B. Tabarrok, M. Farshad, H. Yi,  Finite element formulation of spatially curved and twisted rods,  Computer Methods in Applied Mechanics and Engineering 70 (1988) 275–299.
  * Arnold and Brüls [2007] M. Arnold, O. Brüls,  Convergence of the generalized-$\alpha$ scheme for constrained mechanical systems,  Multibody System Dynamics 18 (2007) 185–202.
  * Brüls and Golinval [2008] O. Brüls, J. C. Golinval,  On the numerical damping of time integrators for coupled mechatronic systems,  Computer Methods in Applied Mechanics and Engineering 197 (2008) 577–588. doi:10.1016/j.cma.2007.08.007.
  * Brüls et al. [2014] O. Brüls, V. Acary, A. Cardona,  Simultaneous enforcement of constraints at position and velocity levels in the nonsmooth generalized-$\alpha$ scheme,  Computer Methods in Applied Mechanics and Engineering 281 (2014) 131–161. doi:10.1016/j.cma.2014.07.025.
  * Bauchau [2010] O. Bauchau, Flexible Multibody Dynamics, Springer Netherlands, 2010.
  * Shabana [2001] A. Shabana, Computational Dynamics, Wiley, 2001.