Assessing Numerical Error in Structural Dynamics Using Energy Balance

Abstract

This work applies the variational principles of Lagrange and Hamilton to the assessment of numerical methods of linear structural analysis. Different numerical methods are used to simulate the behaviour of three structural configurations and benchmarked in their computation of the Lagrangian action integral over time. According to the principle of energy conservation, the difference at each time step between the kinetic and the strain energies must equal the work done by the external forces. By computing this difference, the degree of accuracy of each combination of numerical methods can be assessed. Moreover, it is often difficult to perceive numerical instabilities due to the inherent complexities of the modelled structures. By means of the proposed procedure, these complexities can be globally controlled and visualized in a straightforward way. The paper presents the variational principles to be considered for the collection and computation of the energy-related parameters (kinetic, strain, dissipative, and external work). It then introduces a systematic framework within which the numerical methods can be compared in a qualitative as well as in a quantitative manner. Finally, a series of numerical experiments is conducted using three simple 2D models subjected to the effect of four different dynamic loadings.

1. Introduction

1.1. Targets and Interest of Our Research

Variational mechanics date back as far as the Eighteenth Century, when Leibniz, Euler, Maupertuis, and Lagrange devised the calculus of variations and the principles of least action. This methodology of treating physical phenomena is based on the notion that everything in nature tends to a state of minimal energy [1].

In structural engineering practice, there is a preference to use forces and accelerations rather than energy concepts. Unfortunately, this approach often limits our understanding of the phenomena, as, for example, in the case of earthquakes, damage is a function of the square of the velocity, and not so much of the acceleration [2].

In parallel we will deal with a systematic treatment of the numerical methods which have proliferated since the 1950s with the ever-increasing power of computers. This ceaseless growth in numbers and terminology has given place to a cumbersome mix of mathematics, physics, and computer science that is often difficult to grasp. We have used our work to propose a possible categorization according to the physical qualities which they represent instead of according to their mathematical properties.

1.2. Variational Mechanics

According to the principles of variational mechanics [3], the difference between kinetic energy and strain energy in a structural system equals the applied work due to external forces. In this way, by computing the energy scalars and carefully accounting for this difference at each time step, one should be able to infer the degree of accuracy of a simulation [4].

The correct values should not in any case diverge much from zero, and deviations from this value would give us an idea about how accurate and stable a method is.

1.3. Numerical Methods for Structural Analysis

In previous works published by the authors [5, 6], it was shown how the vast amount of existing numerical methods can be grouped into three main sets according to the kind of physical phenomena they represent and the type of differential equations they discretize: matter integration techniques (partial differential equations), constraint integration techniques (algebraic differential equations), and time integration techniques (ordinary differential equations).

Based on this concept, we have chosen the following matter integration implementations: finite element (FEM), finite differences (FDM), and mass spring systems (MSS). For the constraint integration we will limit ourselves to the constraint reduction (CR) technique, whereas, in the case of time integration we will study the Newmark Beta (NB), Hilber-Hughes-Taylor (HHT), Chung-Hulbert's generalized-alpha (CH), and Wilson Theta (WTH) methods.

All these time integration methods are available in a general-purpose commercial package, so we were able to establish a comparative reference for our own implementations. In the case of matter integration, we implemented our own algorithms from the literature and adapted them to our own purposes, also making a previous benchmark of their results with respect to those obtained by the aforementioned software.

1.4. Numerical Experiments

Three simple structural models under four dynamic loadings will be tested. The influence of the parameters time step, damping ratio, and the number of integration points will be studied.

The work done by the load patterns, along with the internal elastic, kinetic, and dissipative energies, will be computed at each time step and combined together to verify the Hamiltonian energy balance. Its integral through time will provide different values of the total Lagrangian action of the structure-loads system. The deviation from a proposed analytical value, whose computation is straightforward, would account for the level of accuracy of the implementations.

It will be shown how, whether used on single elements or complex systems with more elements, this methodology could be employed as a reference since the value of the action is a simple scalar which is easy to monitor.

2. Variational Mechanics

2.1. Principle of Least Action

In variational mechanics, the Lagrangian functional L, describing the dynamics of a system, is given by

\begin{matrix} L = T - U, \end{matrix}

(1)

where T and U are the kinetic and potential energies of the system, respectively.

According to Hamilton's definition, action S is the integral through the studied time lapse of the Lagrangian

\begin{matrix} S = \int_{t_{1}}^{t_{2}} (T - U) d t . \end{matrix}

(2)

The correct path for a dynamic system is the one for which the value of the action integral is stationary. This leads to a minimization problem which is rooted in the variational principles of Lagrange and Euler.

2.2. Euler-Lagrange Equation and Energy Balance

For a single particle-spring system subject to an external force, the Lagrangian can be written as

\begin{matrix} L (x, \dot{x}) = \frac{1}{2} m \cdot {\dot{x}}^{2} - \frac{1}{2} k \cdot x^{2}, \end{matrix}

(3)

where m is the mass of the particle, k is the stiffness of the spring, x is the instantaneous position, and the superscript dot indicates derivation with respect to time.

From Hamilton's principle of stationary action, and after some variational calculus, the evolution of a physical system is described by the solutions of the forced Euler-Lagrange equation for the action of the system:

\begin{matrix} \frac{d}{d t} \frac{\partial L (x, \dot{x})}{\partial \dot{x}} - \frac{\partial L (x, \dot{x})}{\partial x} = Q_{ext} (x, \dot{x}), \end{matrix}

(4)

where

\frac{\partial L (x, \dot{x})}{\partial x} = \frac{d U (x)}{d x} = \frac{d}{d x} (\frac{1}{2} k \cdot x^{2}),

(5)

\frac{\partial L (x, \dot{x})}{\partial \dot{x}} = - \frac{d T (\dot{x})}{d \dot{x}} = - \frac{d}{d \dot{x}} (\frac{1}{2} m \cdot {\dot{x}}^{2}),

(6)

\begin{matrix} Q_{ext} (x, \dot{x}) = - c \cdot \dot{x} + f_{ext} (t) . \end{matrix}

(7)

Substituting (5), (6), and (7) into (4), and deriving (5) with respect to time, we get the Newtonian classical formulation:

\begin{matrix} \frac{d}{d t} \frac{\partial L}{\partial \dot{x}} - \frac{\partial L}{\partial x} = m \cdot \ddot{x} + c \cdot \dot{x} + k \cdot x - f_{ext} (t) = 0, \end{matrix}

(8)

where the externally applied force f_ext is generally given, and the velocity dependent damping term is a nonconservative force defined in terms of d'Alembert's virtual work [7].

2.3. Kinetic Energy of a System, T

For a structural system under dynamic forces, the above equations are used in a vector-matrix fashion. Each of the points of the structure and its degrees of freedom are represented as terms of a vector, and the mass and stiffness of the whole system are characterized by a matrix. This leads to the following expression for the computation of the kinetic term:

\begin{matrix} T = \frac{1}{2} \cdot {\dot{x}}^{T} \cdot [M] \cdot {\dot{x}}, \end{matrix}

(9)

where ${\dot{x}}$ is the vector of the velocities obtained.

In the present work, the construction of the mass matrix consists of the simple addition of the elements' particular masses in their concurrent nodes (lumped mass matrix).

2.4. Elastic Potential Energy, U

When a body of some material is subject to external forces, its internal structure is deformed. The displacement of these forces in space is the source of performed work.

The scalar value of such work, in order to preserve the balance of energy, must be equal to that of the internal forces in the body (stresses) times the internal displacements within the material (strains).

The energy that is not recoverable is commonly dissipated in the form of heat. However, for the sake of simplicity the scope of this paper will be limited to the elastic range.

In elastic materials, the stored potential strain energy can be accounted for as half of the integral over the volume of the internal strains times the internal stresses. The corresponding formula is [8]

\begin{matrix} U_{el} = \frac{1}{2} \int_{V}^{} {σ}^{T} {ε} d V, \end{matrix}

(10)

where

{σ}^{T} = {\begin{Bmatrix} σ_{x x} & σ_{y y} & σ_{z z} & τ_{x y} & τ_{x z} & τ_{y z} \end{Bmatrix}},

(11)

\begin{matrix} {ε}^{T} = {\begin{Bmatrix} ε_{x x} & ε_{y y} & ε_{z z} & γ_{x y} & γ_{x z} & γ_{y z} \end{Bmatrix}} . \end{matrix}

(12)

In the case of the linearized beam depicted in Figure 1, we can define four kinds of strain energy according to the four main stress components: axial (N), shear (Q), bending moment (M), and torsional moment (T).

Figure 1:

Stress-strain components in a beam. The directions of the infinitesimal strains and stresses are arranged according to the length of the beam.

In Table 1, the formulae for each of these strain energy components are given. The listed expressions can be either a function of the displacements along the beam or of the input forces.

Table 1:

Displacement and force based formulae for the elastic strain energy in a beam.

Displacement		Force

Axial	$U_{A} = \frac{1}{2} \int_{0}^{l} E A \frac{d u^{2}}{d x} d x$	$U_{A} = \frac{1}{2} \int_{0}^{l} \frac{F^{2}}{E A} d x$
Bending	$U_{M} = \frac{1}{2} \int_{0}^{l} E I^{2} d x$	$U_{M} = \frac{1}{2} \int_{0}^{l} \frac{M^{2}}{E I} d x$
Shear	$U_{S} = \frac{1}{2} \int_{0}^{l} A G \frac{d u^{2}}{d x} d x$	$U_{S} = \frac{1}{2} \int_{0}^{l} \frac{F^{2}}{A G} d x$
Torsion	$U_{T} = \frac{1}{2} \int_{0}^{l} G J \frac{d^{2}}{d x} d x$	$U_{T} = \frac{1}{2} \int_{0}^{l} \frac{T^{2}}{G J} d x$

For a structural system, where several elements are combined and attached in n nodes, the equations that establish the behaviour of each node with respect to others are defined in the stiffness matrix [K].

The coefficients that compose this matrix are obtained through the different matter integration methods (FEM, FDM, MSS, BEM, etc.) by solving the above equations in a combination for all three kinds of stresses in all three planes [9].

Finally, in order to compute the total elastic energy U of the system, the following expression:

\begin{matrix} U = \frac{1}{2} {x}^{T} \cdot [K] \cdot {x} \end{matrix}

(13)

is used, where {x} is the vector of the displacements obtained.

2.5. Work Done by Dissipative Forces

In every real structure the existence of damping is a well-known phenomenon whose nature is still not fully understood. In order to incorporate it into a simulation, numerical artefacts are created which can account for the energetic dissipation involved.

In general, a damping matrix [C] is defined which accounts for the dissipative properties of the structural elements. This matrix affects the velocity in Newton's equation (8), behaving as a force which acts against the applied external forces.

The work done by this force can be accounted for by means of the following relation:

\begin{matrix} R = \frac{1}{2} {x}^{T} \cdot [C] \cdot {\dot{x}} . \end{matrix}

(14)

The simplest model for dissipation in structural dynamics is due to Lord Rayleigh and is known as “linear damping,” “Rayleigh damping,” or “classical damping.” In this idealization, the damping matrix is assumed to be a linear combination of the stiffness and the mass matrices. Despite the numerous criticisms this model receives, it is still widely used for its convenience when combined with the modal analysis procedure [10]. Once the stiffness matrix [K] and the mass matrix [M] are ready, the damping matrix [C] can be defined as follows:

\begin{matrix} [C] = α \cdot [M] + β \cdot [K] . \end{matrix}

(15)

The value of the mass and stiffness coefficients is determined by the solution of the eigenvalues of the [K] matrix [10].

2.6. Work Done by External Forces

The total work exerted over a structure by the external applied forces can also be represented in a vectorial fashion as

\begin{matrix} W_{ext} = \frac{1}{2} {F_{ext}} \cdot {x} . \end{matrix}

(16)

In general, the vector of the external forces is known for each time step.

2.7. Total Action of the System, Energy Balance, and the Lagrange-d'Alembert Principle

In order to account for the correctness of a simulation, we can utilize the Lagrange-d'Alembert principle [11], which establishes the following relation:

\begin{matrix} δ \int_{t_{1}}^{t_{2}} L d t + \int_{t_{1}}^{t_{2}} F_{ext} δ x d t = 0 . \end{matrix}

(17)

If we withdraw the variation operator and rearrange the terms, this leads to

\begin{matrix} \int_{t_{1}}^{t_{2}} L d t = - \int_{t_{1}}^{t_{2}} F_{ext} x d t, \end{matrix}

(18)

which, in discrete form, yields

\begin{matrix} \sum_{t_{1}}^{t_{2}} L d t = - \sum_{t_{1}}^{t_{2}} F_{ext} x d t . \end{matrix}

(19)

Having defined previously each of the terms, we can now write an elementary formula from which the degree of exactness of a simulation can be estimated:

\begin{matrix} \sum_{t_{1}}^{t_{2}} (T - U) d t = - \sum_{t_{1}}^{t_{2}} (W_{ext}) d t . \end{matrix}

(20)

This is basically the computation of an energy balance, where the Hamiltonian action is treated, in its discrete form, as an average over time of each instantaneous Lagrangian. In order to account for the external forces involved, we also integrate their work over time. According to d'Alembert's principle, these two measures should be equal when internal dissipative forces (hysteretic damping) are not present.

It is noteworthy how (20) also represents the first law of thermodynamics, where the value of heat (in this case added to the structure) can be associated with the kinetic energy of the system.

Any divergence from this equality gives a measure as to how inaccurate a numerical method is by means of a single value, without the need to find simplified analytical models whose assumptions rarely fit the real problems of engineering practice.

3. Numerical Methods

As explained in previous works by the authors [5, 6], for the simulation of structural dynamics three different physical notions need to be integrated: time, matter, and kinematic constraints.

Figure 2 shows a schematic of the process involved and its relation to the physical concepts of time, matter and constraints. Also for convenience, a list of abbreviations for existing numerical methods is provided in Table 2.

Table 2:

Abbreviations employed for some numerical methods.

Abbreviation	Method

FEM	Finite element method
FDM	Finite differences method
FVM	Finite volumes method
MSS	Mass spring systems
SPH	Smoothed particle hydrodynamics
PU	Partition of unity
MLS	Mean Least squares
PM	Penalty method
LM	Lagrange multipliers
GC	Generalized coordinates
UK	Udwadia-Kalaba
IB	Impulse based
RK	Runge-Kutta
LAM	Linear acceleration method
LF	Leapfrog
IE	Implicit Euler
BDF	Backward difference formula
NB	Newmark-Beta
HHT	Hilber-Hughes-Taylor

Figure 2:

Schematic of some numerical methods and their associated physical notions. The arrow represents a possible sequence of methods for a dynamics simulation.

3.1. Matter Integration

To describe the dynamics of matter we have an infinite number of degrees of freedom because the particles that compose it can have arbitrary displacements with respect to each other. Such systems are described using partial differential equations (PDEs), where the time and spatial coordinates are related. These general partial differential equations, which are applicable to any solid or fluid material, are derived from the constitutive laws which apply to the material.

For the solution of these equations, two different approaches can be taken in order to control the number of degrees of freedom (i.e., discretize): creating a mesh where the material displacements are limited (mesh based methods) or establishing the equations in the form of potential functions so that they compose a system of particles that regulate each other (mesh free methods) [12, 13].

We have limited our study to three mesh based methods with different discretization schemes: finite element method (FEM), finite differences method (FDM), and a mass spring system (MSS).

For the general computation of nodal displacements and rotations, a framework employing the direct stiffness method (DSM) was prepared [14].

In our case, where beam elements were used, the analytical solution of Bernoulli-Euler is lumped into local element matrices that are ultimately assembled into a global stiffness matrix [8].

For FEM implementation, the description of the elastic deformation of the beam is based on a Hermite interpolation polynomial obtained from [3].

FDM establishes the relations between stations along the beam as a sequence of equations that form a linear system which is easily invertible [9, 15].

MSS is somewhat more complex since it requires previous discretization of the beam into a set of connected tetrahedra, but from the point of view of physics its results are clearer as the assumptions are that the nodes are simply connected by bars with a characteristic Young's modulus and area [16]. Some adjustments had to be made to the position of the masses in the cross section so that the inertia of the section would match the value assigned in the polynomial-based methods.

The global nodal displacements and rotations computed by means of DSM were transformed ultimately into local coordinates and served as input variables for each of the three methods above.

3.2. Kinematic Constraints Integration

In order to obtain the internal strains and stresses, bodies are subject to kinematic constraints. In this way a set of differential algebraic equations (DAEs), comprising the PDEs mentioned in the previous chapter, is assembled in the respective stiffness, mass, and damping matrices:

\begin{matrix} K g = [\begin{bmatrix} \cdot & \notin & \cdot & \cdot & \cdot & \cdot \\ \notin & \notin & \notin & \notin & \notin & \notin \\ \cdot & \notin & \cdot & \cdot & \cdot & \cdot \\ \cdot & \notin & \cdot & \cdot & \cdot & \cdot \\ \cdot & \notin & \cdot & \cdot & \cdot & \cdot \end{bmatrix}] ⟹ K g_{redd} = [\begin{bmatrix} \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \end{bmatrix}] . \end{matrix}

(21)

Equation (21) shows DAE integration by constraint reduction (CR). The global stiffness matrix of the system is made nonsingular by symmetrically subtracting the columns and rows corresponding to the constrained degrees of freedom (Kg: original matrix; Kg_redd: modified reduced nonsingular matrix). Consider

\begin{matrix} K g = [\begin{bmatrix} \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{bmatrix}] ⟹ K g_{ext} = [\begin{bmatrix} \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot \end{bmatrix}] [\begin{bmatrix} 1 & \cdot \\ \cdot & 1 \\ \cdot & \cdot \\ \cdot & \cdot \\ \cdot & \cdot \end{bmatrix}] [\begin{bmatrix} 1 & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & 1 & \cdot & \cdot & \cdot & \cdot \end{bmatrix}] [\begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix}] . \end{matrix}

(22)

Equation (22) shows DAE integration by Lagrange multipliers (LM). The global stiffness matrix is made nonsingular by symmetrically adding columns and rows where unit values are placed in the location of the constrained degrees of freedom (Kg: original matrix; Kg_ext: modified extended nonsingular matrix). Consider

\begin{matrix} K g = [\begin{bmatrix} \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{bmatrix}] ⟹ K g_{sc} = [\begin{bmatrix} \infty & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \infty & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \infty \end{bmatrix}] . \end{matrix}

(23)

Equation (23) shows DAE integration by penalty method (PM). The singularity of the global stiffness matrix is treated by scaling the diagonal elements of the constrained degrees of freedom with a very large number (Kg: orginal matrix; Kg_sc: modified scaled nonsingular matrix).

There is vast a amount of literature dedicated to explain different schema from which stable, accurate, and faster formulations can be constructed [12]. In general terms, Newton's equation is manipulated in one way or another to isolate one of its terms. The possibilities are to do it either at the acceleration level, at the velocity level, or at the position level.

In the most common acceleration level schemes, predominance is given to the methods of constraint reduction (CR), Lagrange multipliers (LM), and penalty method (PM).

In all three methods the strategy is to alter the stiffness, damping, and mass matrices in such a way that they become invertible (after assembly, the stiffness matrix is symmetrical and singular).

By means of CR one simply removes from the stiffness and mass matrices the rows and columns corresponding to the constrained degrees of freedom, hence reducing it in size.

Similarly, the LM method increases the size of the system of equations in as many rows and columns as there are constrained degrees of freedom in the model. These “attached” vectors are zero in all their elements, except for those corresponding to the aforementioned DOFs, where a value in the range [0,1) is inserted.

Regarding the Penalty Method, the corresponding elements of the diagonal are replaced with values that have several orders of magnitude larger than those initially assembled.

Equation (23) and Figures 3 and 4 provide a visualization of these methodologies as they are commonly implemented. Further analysis on their advantages and drawbacks can be found in [17].

Figure 3:

Geometry of the three models. Dimensions in cm.

Figure 4:

Frequency response functions for the three tested models.

3.3. Time Integration

Integration over time in a structural dynamics simulation reduces to the solution for each time step of an ordinary differential equation (ODE). The first possible classification for ODEs solvers distinguishes between explicit, implicit, and hybrid methods. This division arises as a consequence of the so-called numerical stiffness. This phenomenon forces the size of the adopted time step to be so small that convergence is never reached; or otherwise, the adopted time steps are so large that the simulation becomes unstable. This can be caused by the physical characteristics of the system (components with large differences in their masses, stiffness, and/or damping). However, in many other instances, stiffness is numerically induced due to either the discretization process, the large number of components and equations of motion, or sudden or accumulated violations in the constraint conditions.

Explicit methods present this kind of problem. The advantage of implicit methods is that they are usually more stable for solving a stiff equation, meaning also that a larger step size can be used. Nevertheless, extra computations need to be performed internally, thus requiring extra time.

For our comparison we have used the methods available in a general purpose commercial package: Newmark Beta (NB), Wilson Theta (WTH), Hilbert Hughes Taylor (HHT), and Chung and Hulbert (CH). We implemented our algorithms from [17 –21]. The results obtained were in very good agreement with those of the commercial package.

4. Numerical Experiments and Results

In this chapter we provide the results of our numerical experiments, where several combinations of methods were used in diverse simulations. Three different specimens of increasing complexity were tested, and some engineering relevant parameters affecting each numerical method were systematically studied (time step, damping ratio, and number of integration points).

In order to avoid excessive complexity, the specimens were treated as 2D models and kept within the elastic range, considering the shear effects in deformation to be negligible.

4.1. The Studied Specimens

As mentioned above, and for the sake of simplicity, we omitted material and geometrical nonlinearities from our analyses.

The material and geometric properties shown in Table 3 are common in engineering practice, with values similar to those of a 200 × 200 × 2 mm hollow extruded steel bar.

Table 3:

Properties of the beam elements composing the specimens.

Parameter
Area, A	144 cm²
Modulus of inertia, I	7872 cm⁴
Modulus of elasticity, E	21000 KN/cm²
Shear modulus, G	8076,92 KN/cm²
Density, d	7.892E – 8 KN/cm³

The geometric configuration of each model is displayed in Figure 3, in order of increasing complexity.

Notwithstanding the obvious resemblance to a typical structural engineering application, the generality of our work and its potential use in the simulation of any mechanical object regardless of size or shape should be noted. However, the choice of parameters to characterize the strain energy must be judicious in order to achieve adequate stiffness and mass matrices.

Model A. The simplest model of choice for our research was a 387,5 cm long cantilever column under a lateral loading acting in its tip. The cantilever model is extensively utilized for the validation of numerical methods in the literature. It is composed of two elements, each half the length of the column.

Model B. A natural extension to this model from the structural engineering point of view is a simple moment frame, with identical geometrical and mechanical properties for each beam element as in the previous case. The load F is applied to the upper left corner.

Model C. The more complex three bay-four storey frame is also shown in Figure 3. Its properties are again those presented in Table 3, and load F is also applied to the upper left corner.

In order to represent the structural dissipative behaviour, Rayleigh damping was implemented according to [10]. It is based on modal analysis and uses the first two natural frequencies of the structure under study. The ones applicable to our models are listed in Table 4.

Table 4:

Modal frequencies for damping characterization.

Model	Mode	Frequency (Hz)

A	1st	12.79
A	2nd	64.44
B	1st	11.37
B	2nd	33.52
C	1st	2.71
C	2nd	8.69

For comparison purposes, a frequency response function was computed for all three models. Its values are in agreement with those of the modal analysis of Table 4 as can be seen in Figure 4. It can be inferred from this figure that the more complex model C has the highest sensitivity to low frequencies, whereas models A and B should behave similarly as they have their strongest response to similar frequency values.

4.2. Transient Input Forces

A load F of 10 KN applied to the tip of each specimen was scaled at each time step with an input signal of variable amplitude.

As presented in Figures 5, 6, 7, and 8, four input signals were devised in order to stimulate the loading of our system: a simple sine function, a simple sine function suddenly interrupted, an incremental triangular function, and a ramp pulse, all of them spanning through five seconds.

Figure 5:

Sine function, two cycles: f = 0.4 Hz, T = 2.5 s.

Figure 6:

Sine function, one cycle, then free vibration: f = 0,4 Hz, T = 2,5 s.

Figure 7:

Incremental triangular function: f = 1.2 Hz, T = 0.83 s.

Figure 8:

Ramp pulse: F = 0.625 Hz, T = 1.6 s.

A sine function with such a low frequency is seldom encountered in engineering practice, but allows for the calibration and tuning of the combined methods given its smoothness and clarity.

For the second signal, after completion of the first period it is interrupted abruptly in order to allow for free vibration of the system. The point of interruption, in zero amplitude, allows for observation of the effect of kinetic energy on the simulation.

The incremental triangular function was constructed in order to account for earthquake engineering regulations, where sudden changes and peaks are to be simulated.

Regarding the last pulse, it makes it possible to compare the performance of the numerical methods simulating free vibration and the effect of resonance.

4.3. Parametric Sensitivity Study

The significant parameters involved in the numerical computations have been iteratively modified in order to assess their influence in the simulations. For each type of integration the following parameters were studied:

Time integration:

time step influence,

damping ratio influence.

Matter integration:

number of integration points along the beam element.

Constraint integration:

no comparison was available, as only the constraint reduction technique is implemented in the reference software.

Table 5 shows the values used for the characteristic parameters of each numerical method in all the simulations.

Table 5:

Time integration parameters.

Method	Alpha	Beta	Gamma	Alpha-m	Theta

NB	—	1/4	1/2	—	—
WTH	—	—	—	—	1.4
HHT	−1/3	0.444	0.8333	—	—
CH	−1/3	1/4	1/2	−1/10	—

These values were not the subject of our study and were fixed according to recommended values from the literature [18 –21]. It is important to note that Chung-Hulbert's method (also known as Generalized-Alpha) under certain combinations of parameters includes previous ones, whose performances are, according to [21], less accurate when low frequency excitation is present.

4.4. Methodology: Energy Computation of a Simulation

The evaluation of instantaneous energetic magnitudes provides a very holistic hindsight into the behaviour of a simulation, which is qualitatively superior to that of the displacement domain to which time history analysis is traditionally limited.

Besides, in the case of the single cantilever beam choosing the tip as the observed target is generally straightforward, but for more complex arrangements like, say, models B and C, this is not so trivial. The common choice of a “representative point” (the centre of mass of each storey, conversion to SDOF, etc.) has a definition which is always difficult and elaborate.

As an example, simple observation of the displacement behaviour of the tip of Model A would mislead the analyst to the conclusion that the results for signal 4 in Figure 9 are better approximations than those for signal 3, as the displacement values seem to be closer to the analytical ones.

Figure 9:

Model A. Time history analysis of the displacement of the tip. Chung-Hulbert method, generalized α = – 0.1, dt = 0.0025, and damping ratio = 1%. Analysis in the displacement domain is limited to the observation of single nodes of the system and can be misleading, as smooth displacement behaviour often conceals energy conservation violations.

Nevertheless, this can be proven to be less accurate than expected. Figure 10 shows the same simulation in the energy domain, computing some operators of the different terms from Section 2. The application of (20) appears in Figure 10 as W + T – U – R. From this operator one can obtain that, on average, under signal 4 the simulation “creates” +0.69 N cm of spurious energy on each time step, whereas the same model, under signal 3—visibly more flurry in the displacement domain, “absorbs” −1.12 N cm from nowhere. In terms of absolute value, the first is closer to zero, apparently still showing a better approximation for signal 4. However, a rigorous computation should also take into account that the total amount of work applied by signal 3 is, on average, three times larger than that of signal 4. It is not equivalent to a large average deviation from zero with large values as it is with smaller ones. The formulation of an independent normalization parameter will be provided.

Figure 10:

Model A. Time history for the variation of different energy operators. Chung-Hulbert method, generalized α = – 0.1, dt = 0.0025, and damping ratio = 1%. Analysis in the energy domain not only gives a global view of the structural behaviour but also accounts for possible points of numerical instability and irregularities of the simulation.

To define our measure of error we use (9), (13), (14), and (16) at each time step to compute the respective instantaneous values of kinetic energy (T), strain energy (U), dissipative energy (R), and external work (W).

Our methodology, based on (20), uses the Hamiltonian action integral minus the average over time of the work due to the externally applied forces, thus measuring the difference to zero.

Moving the Lagrangian (T-U) term to the right hand side and including the damping term, we obtain:

\begin{matrix} ϵ (t) = W_{ext} (t) + T (t) - U (t) - R (t), \end{matrix}

(24)

whose discrete integral in time gives

\begin{matrix} \sum_{t_{1}}^{t_{2}} ϵ (t) d t = \sum_{t_{1}}^{t_{2}} [W_{ext} (t) + (T (t) - U (t) - R (t))] d t . \end{matrix}

(25)

However, this value is, by itself, not yet very representative. As mentioned before, in the case of large values of external work, epsilon also yields large values. For this reason, a normalization parameter was devised. It is based on the total work done by the external forces, but computed independently from the displacements and based on (13) and (16). It is obtained by isolating the displacement vector in each equation:

{x}^{2} = 2 U \cdot {[K]}^{- 1},

(26)

\begin{matrix} {x} = 2 W_{ext} \cdot {F_{ext}}^{- 1} . \end{matrix}

(27)

By combining (26) and (27) it is now possible to define our reference input work parameter:

\begin{matrix} W_{ref} = 2 \frac{W_{ext}^{2}}{U} = {F_{ext}}^{T} {⌊ K ⌋}^{- 1} \cdot {F_{ext}} . \end{matrix}

(28)

This reference work is completely self-contained and does not rely on the numerical method used to do the simulation, as the vector of the external forces and the stiffness matrix are given data.

Once this value is defined, it is possible to establish a measure of the error that is consistent regardless of the numerical method, applied force, or desired model of study.

The applied formula for each of our comparisons is then as follows:

\begin{matrix} Error = \frac{100 \cdot \sum_{t_{0}}^{t_{1}} ϵ d t}{\sum_{t_{0}}^{t_{1}} W_{ref} d t} . \end{matrix}

(29)

Other options for the value of epsilon are also available. Similarly, one could compute (20) using the Hamiltonian (T + U + R) and subtracting it from the applied work. Its time history is shown in Figure 10 as W – (T + U + R). This operator provides a lower bound for the evolution of the Lagrangian (most clearly visible for signal 4), as it balances the kinetic energy of the system against the potential and the dissipative energies. Its evolution in time gives information about whether the absolute value of the kinetic term is overestimated at each step. Given that the mass is kept constant this operator permits to verify that instantaneous velocities are computed correctly.

Another possible option for the value of epsilon is to calculate the instantaneous increment of the Hamiltonian, d(T + U + R). In systems where the energy is constant, this value should be zero, but it is rarely the case in practical applications. Its main interest resides in the detection of segments in the simulation where the smooth transition from one time step to the next is lost.

One could also define the epsilon on each time step as the difference between the time-dependent calculated work and our presented analytical reference work (W_ref – W). In a way, this computation appears the most precise, as the involved terms are of the same kind and the reference work is derived from a numerically neutral relationship. Apart from the possible error in the inversion of the stiffness matrix, the term W_ref is immune to the fluctuations caused by the time integrators. Still, this operator is not fully satisfactory. As the possible errors in the instantaneous work only depend on the computed displacement, its time history only provides information about irregularities in this matter.

The choice, then, of the Lagrangian (minus the damping energy when applicable) to balance the external work seems the most appropriate. Not only is its time history a valuable source of information for the analysis of irregularities in a simulation but also its integral in time provides a single scalar whose value should be zero. Given that the energetic terms are all positive, a positive value of this integral can only be caused by an average overestimation of the kinetic term (i.e., the velocities) against the displacements. Similarly, a negative value tells us to what degree displacements are unbalanced against velocities, as the internal potential energy is a direct function of them.

4.5. Numerical Results: Influence of Time Step

According to the methodology exposed above, a thorough parametric study was carried using the three models of choice. Figures 11, 12, and 13 present the values obtained from iteratively modifying the time step between values of 0.00125 s and 0.15 s for each numerical method, with a constant damping ratio value of 2%.

Figure 11:

Energy error analysis. Model A: influence of time step size. Damping ratio = 2%.

Figure 12:

Energy error analysis. Model B: influence of time step size. Damping ratio = 2%.

Figure 13:

Energy error analysis. Model C: influence of time step size. Damping ratio = 2%.

As opposed to the analyst's intuition, in spite of dealing with linear models we obtained curves that vary significantly from one method to another. Nevertheless, and as expected, this divergence is more pronounced with larger time steps and also increases with the complexity of the model.

The character positive or negative of the value of the error also provides a valuable source of information, as it tells us when the internal strain energy is larger or smaller than the sum of the kinetic energy plus the external work. As this term is dependent on velocity, it shows when the kinetic term is overestimated or underestimated. In other words, the higher the decoupling between velocity and displacement, the further the simulation is from correctness.

When the time step is larger, it affects the velocity, which loses or gains in phase with the normal modes of the structure and with the input signal. In these cases the simulation might either dissipate or absorb energy artificially. This explains the ripple around the abscissa presented by all the methods in all the simulations.

In terms of evaluation of the particular methods, it is commonly accepted that CH has better performance than the others, as it gives the analyst control over the numerical damping for high frequencies without loss of accuracy. As the sensitivity to those parameters was not within the scope of this study, we cannot give a view about such effect, but we can point out how, in general, in this configuration they all show fairly similar results, only diverging significantly for larger and impractical time steps. Although all of them are of the implicit type, meaning unconditional stability regardless of the time step size, our results show how this set of methods in general tends to sacrifice energy conservation. In most linear structural dynamics problems it is still not an issue, but for the analysis of nonlinear situations we strongly recommend the use of more modern integrators of the symplectic type, as those described in [22].

4.6. Numerical Results: Influence of the Damping Ratio

It should be noted that the damping considered in our experiments is of an external nature, given the fact that no material nonlinearities have been taken in to consideration. The corresponding Rayleigh mass and stiffness coefficients defined in (27) were obtained according to [10]. Figure 14 shows the relation of these values to the models used in the study.

Figure 14:

Rayleigh damping coefficients. Their values are directly proportional to the value of the chosen damping ratio. For higher frequencies of the model, the value of the mass coefficient is higher, and vice versa for the stiffness coefficient.

The sensitivity of the numerical methods to variations in the damping ratio is presented in Figures 15, 16, and 17. For all three models the range of study was fixed between 0% to 10% of critical damping. In general, this is sufficient for all the methods to reach their asymptotic limit in almost every simulation. For models A and B a value of 2% of damping suffices to achieve stable behaviour with an error of less than 0.3%, which can be considered very acceptable.

Figure 15:

Energy error analysis. Model A: influence of damping ratio. Time step = 0.01 s.

Figure 16:

Energy error analysis. Model B: influence ofdamping ratio. Time step = 0.01 s.

Figure 17:

Energy error analysis. Model C: influence ofdamping ratio. Time step = 0.01 s.

4.7. Numerical Results: Influence of the Number of Integration Points for Matter Integration Methods

For the study of the matter integration techniques a similar approach based on the variational principle of action was adopted. However, here the definition of a reference parameter W_ref was not required as the analytical solution for beam elements is available applying the concepts of Section 2.4.

Instead, we computed a global action term, whose units are also those of angular momentum. It is defined as

\begin{matrix} S = \int_{t}^{} U_{elas} (t) d t . \end{matrix}

(30)

The charts in Figure 18 were made by computing the instantaneous value of internal work corresponding to each different numerical method and averaging it over time. The applied transient force was signal 1.

Figure 18:

Comparison of error in angular momentum computation for matter integration methods against the number of integration points. Analytical (ANA) versus finite differences (FDM) versus finite element (FEM) versus mass spring system (MSS).

The measure of the error was computed as a percentage of the difference to the analytical value. A positive error indicates numerical spurious dissipation of energy, whereas a negative error stands for artificial energy creation.

As expected, for an increasing number of integration points the methods converge towards the analytical value. However, they do it in an asymptotic fashion, reaching an almost flat parallel value after about 25 integration points. In general, the obtained error values remain below 5% for all the methods, which is completely acceptable in practice.

Interestingly, FEM presents the best behaviour only for the simple cantilever beam, creating spurious strain energy for the other two models.

FDM and MSS tend to dissipate energy in all cases, which means that, in general, they result in an underestimated value of displacement by about 2%, remaining on the unsafe side.

5. Discussion and Future Work

The total Hamiltonian actions of three systems under transient loadings have been computed for each possible combination of methods. A comparison was made on the basis of energy principles. It was shown how variational principles and an energetic norm can be employed in an easy and efficient manner to benchmark and assess the accuracy and stability of different implementations. The scheme provided, tested in three simple examples, is trivially extensible to more complex systems where more elements are present. The advantage of this approach is that it allows for the monitoring of the global behaviour by means of one simple scalar. Also, a conceptual framework for the classification and treatment of numerical methods, grouping them into time, matter, and constraint integrators, was used for systematic analysis of the results. It was also shown how methods of a different nature and concept can be compared using the same theoretical background. The accuracy and good performance of time and matter integration methods are generally taken for granted, as it is difficult, in the displacement domain, to assess them with certainty. Future work will aim at the application of the same methodology in material nonlinear analysis and in structures of higher complexity.

References

Lanczos

, The Variational Principles of Mechanics, University of Toronto Press, Toronto, Canada, 1970.

Housner

G. W.

, “Limit design of structures to resist earthquakes,” in Proceedings of the World Conference on Earthquake Engineering, Earthquake Engineering Research Institute, 1956.

Wunderlich

and Pilkey

, Mechanics of Structures. Variational and Computational Methods, CRC Press, New York, NY, USA, 2002.

Bugeda

, Estimación y Corrección Del Error en el Análisis Estructural Por el MEF, Monografia 9, CIMNE, Barcelona, Spain, 1991.

Andújar

Roset

, and Kilar

, “Beyond FEM: overview on physics simulation tools for structural engineers,” Technics Technologies Education Management, vol. 6, no. 3, pp. 555–571, 2011.

Andujar

Roset

, and Kilar

, “Interdisciplinary approach to numerical methods for structural dynamics,” World Applied Sciences Journal, vol. 14, no. 8, pp. 1046–1053, 2011.

Wang

Q. A.

and Wang

, “Is it possible to formulate least action principle for dissipative systems?” Cornell University Library, 2012, http://arXiv.org/abs/1201.6309.

Argyris

J. H.

, Energy Theorems and Structural Analysis, Butterworths Scientific Publications, London, UK, 1960.

Aguiar

, Analisis Estatico de Estructuras, Colegio de Ingenieros Civiles de Pichincha, Quito, Ecuador, 2005.

10.

Adhikari

, Damping models for structural vibration [Ph.D. thesis], Trinity College, Cambridge, UK, 2000.

11.

Bullo

and Lewis

A. D.

, Geometric Control of Mechanical Systems, Texts in Applied Mathematics, Springer, New York, NY, USA, 2004.

12.

Nealen

Müller

Keiser

Boxerman

, and Carlson

, “Physically based deformable models in computer graphics,” Computer Graphics Forum, vol. 25, no. 4, pp. 809–836, 2006.

13.

Liu

G. R.

, Mesh Free Methods: Moving Beyond Finite Element Methods, CRC Press, New York, NY, USA, 2003.

14.

Wasfy

T. M.

and Noor

A. K.

, “Computational strategies for flexible multibody systems,” Applied Mechanics Reviews, vol. 56, no. 6, pp. 553–613, 2003.

15.

Przemieniecki

J. S.

, Theory of Matrix Structural Analysis, McGraw-Hill, New York, NY, USA, 1968.

16.

Smit

G. D.

, Numerical Solution of Partial Differential Equations by Finite Difference Methods, Oxford Applied Mathematics and Computing Science Series, Oxford University Press, Oxford, UK, 2nd edition, 1978.

17.

Wilson

E. L.

, Static and Dynamic Analysis of Structures, Computers and Structures, New York, NY, USA, 2004.

18.

Bathe

K. J.

, Finite Element Procedures in Engineering Analysis, Prentice Hall, New York, NY, USA, 1995.

19.

Newmark

N. M.

, “A method of computation for structural dynamics,” ASCE Journal of Engineering Mechanics Division, vol. 85, no. 3, pp. 67–94, 1959.

20.

Hilber

H. M.

Hughes

T. J. R.

, and Taylor

R. L.

, “Improved numerical dissipation for the time integration algorithms in structural dynamics,” Earthquake Engineering and Structural Dynamics, vol. 5, pp. 282–292, 1977.

21.

Chung

and Hulbert

G. M.

, “A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method,” ASME Journal of Applied Mechanics, vol. 60, no. 2, pp. 371–375, 1993.

22.

Kuhl

and Crisfield

M. A.

, “Energy-conserving and decaying algorithms in non-linear structural dynamics,” International Journal for Numerical Methods in Engineering, vol. 45, no. 5, pp. 569–599, 1999.