On the Nonlinear Transient Analysis of Structures

Abstract

This chapter consists of four sub-chapters the first of them is the introduction. The second sub-chapter discusses the suitability of the explicit and implicit methods for seismic analysis of buildings exposed to a substantially nonlinear response. Taking into account a duration typical for earthquakes from several seconds to about twenty seconds, it seems that the implicit method is more suitable. Because of the shape of accelerograms, it is needed to use quite short time steps. Both methods were compared for a heavily nonlinear response, typical for seismicity. The third subchapter discusses the analysis of the time-dependent behavior of concrete using nonlinear time-dependent material models during fast dynamic actions. Huge plastic yielding or damage of the material is reached. The fourth sub-chapter discusses sources of dynamical damping that are needed in earthquake analysis. This subchapter aims to take into account the influence of different kinds of damping in the dynamic analysis of structures. More advanced analysis taking into account the various sources of damping, i.e., material viscosity, plasticity/damage, and friction in connections and supports is strongly recommended in the earthquake analysis.

1. Introduction

In the last decades, a big development in the field of nonlinear and time-dependent material models, including those that are suitable for seismic analysis, were achieved. But to take advantage of these developments, the use of spectral analysis must be abandoned and a method of direct integration of the equation of motion must be used (also called the transient or time history analysis). Several aspects of this area, including the comparison of the explicit and implicit method in seismic analysis [1, 2], suitable material models, and sources of dynamic damping are discussed in this chapter. Subchapter 2 deal with several aspects of the transient dynamic approach in seismic analysis. Subchapter 3 compared the response to seismic loading of various material models, specifically the Drucker-Prager plasticity model [3] and Mazars damage model [4]. Subchapter 4 showed different types of dynamic damping in structural analysis which should not be neglected in transient analyses.

2. Numerical methods of direct integration of equation of motion and their application in seismic analysis

This subchapter aims to compare the suitability of the explicit and implicit methods in the nonlinear earthquake analysis from the point of view of computational time and accuracy. We implemented these methods in RFEM [5] software and we used Ansys software [6] to compare and validate the results. In the following text let us introduce the numerical methods used in the study. The equation of motion of a discrete model of a structure subjected to a dynamic load can be written

$$\mathbf{M}\ddot{\mathbf{u}}\left(t\right)+\mathbf{C}\dot{\mathbf{u}}\left(t\right)+{\mathbf{f}}^{\mathrm{int}}\left(t\right)=\mathbf{Ma}\left(t\right)+\mathbf{Cv}\left(t\right)+{\mathbf{f}}^{\mathrm{int}}\left(t\right)={\mathbf{f}}^{\mathit{ext}}\left(t\right)$$E1

where $\mathbf{M},\mathbf{C}$ are matrices of mass and damping respectively, ${\mathbf{f}}^{\mathit{ext}}$ and ${\mathbf{f}}^{\mathrm{int}}$ are the vectors of the external and internal forces respectively and $\mathbf{u},\mathbf{v}=\dot{\mathbf{u}},\mathbf{a}=\ddot{\mathbf{u}}$ are displacement, velocity, and acceleration, respectively.

The numerical process of direct integration of the equation of motion (1) in a time interval $t\in ⟨t_0,t_m⟩$ is called transient analysis. Time is discretized to finite number of time instants $t_0,t_1,\dots ,t_m$. The distance between individual time instants $\mathrm{\Delta }{t}_i=t_i-t_{i-1}$ is called the time step. The lengths of time steps $\mathrm{\Delta }{t}_i$ influences the accuracy, stability and speed of the numerical solution. At the time $t=0$ an initial conditions $\mathbf{u}\left(t_0\right)={\mathbf{u}}_0$, $\mathbf{v}\left(t_0\right)={\mathbf{v}}_0$ have to be defined. Eq. (1) then can be read as follows:

$$\mathbf{M}{\ddot{\mathbf{u}}}_i+\mathbf{C}{\dot{\mathbf{u}}}_i+{\mathbf{f}}_i^{\mathrm{int}}={\mathbf{Ma}}_i+{\mathbf{Cv}}_i+{\mathbf{f}}_i^{\mathrm{int}}={\mathbf{f}}_i^{\mathit{ext}}$$E2

2.1 Explicit method: Central differences

In explicit methods, we introduce explicit assumption about the course of displacement in interval $⟨t_i,t_{i+1}⟩$ and the knowledge of the displacement vector ${\mathbf{u}}_i$ at time instant $t_i$, and we calculate vectors ${\mathbf{u}}_{i+1},{\mathbf{v}}_{i+1},{\mathbf{a}}_{i+1}$ from (2).

The numerical integration of differential equations uses the substitution of the derivative of the displacement with respect to time. in (2) by

$${\mathbf{v}}_i=\frac{1}{2\mathrm{\Delta }{t}_i}\left({\mathbf{u}}_{i+1}-{\mathbf{u}}_{i-1}\right)$$E3

$${\mathbf{a}}_i=\frac{1}{\mathrm{\Delta }{t_i}^2}\left({\mathbf{u}}_{i+1}-2{\mathbf{u}}_i+{\mathbf{u}}_{i-1}\right)$$E4

we obtain a recurrent formula for ${\mathbf{u}}_i$

$$\left(\frac{1}{\mathrm{\Delta }{t_i}^2}\mathbf{M}+\frac{1}{2\mathrm{\Delta }{t}_i}\mathbf{C}\right){\mathbf{u}}_{i+1}={\mathbf{f}}_i^{\mathit{ext}}-{\mathbf{f}}_i^{\mathrm{int}}+\frac{2}{\mathrm{\Delta }{t_i}^2}{\mathbf{Mu}}_i-\left(\frac{1}{\mathrm{\Delta }{t_i}^2}\mathbf{M}-\frac{1}{2\mathrm{\Delta }{t}_i}\mathbf{C}\right){\mathbf{u}}_{i-1}$$E5

The method has all the advantages of explicit methods as long as $\left[\mathbf{C}\right]=\left[\mathbf{0}\right]$ or $\left[\mathbf{C}\right]=\alpha \left[\mathbf{M}\right]$, where α is the mass coefficient of the Rayleigh damping. The most effectivity is achieved when the mass matrix $\mathbf{M}$ is diagonal. The explicit methods are conditionally stable. The time steps must satisfy the following condition:

$$\mathrm{\Delta }{t}_i\le \frac{T_n}{\pi }$$E6

where $T_n$ is the smallest vibration period.

2.2 Implicit method: Newmark

Implicit methods in transient dynamics are based on the equation of motion (2) at the time instant $t_i$. The numerical solution of the system is performed step by step by the following text.

$${\mathbf{a}}_i=f\left(t_i,{\mathbf{f}}_i,{\mathbf{u}}_{i-1},{\mathbf{v}}_{i-1},{\mathbf{a}}_{i-1}\right)$$E7

with the necessity to evaluate the acceleration at the start of the motion at time $t_0$ from the equation of motion (2)

$${\mathbf{Ma}}_0+{\mathbf{Cv}}_0+{\mathbf{Ku}}_0={\mathbf{f}}_0$$E8

The basic formula of the Newmark method that specifies the relations between displacement, velocity, and acceleration vectors have the following form:

$${\mathbf{u}}_i={\mathbf{u}}_{i-1}+\mathrm{\Delta }{t}_i{\mathbf{v}}_{i-1}+\left(\frac{1}{2}-\beta \right)\mathrm{\Delta }{t_i}^2{\mathbf{a}}_{i-1}+\beta \mathrm{\Delta }{t_i}^2{\mathbf{a}}_i$$E9

$${\mathbf{v}}_i={\mathbf{v}}_{i-1}+\left(1-\gamma \right)\mathrm{\Delta }{t}_i{\mathbf{a}}_{i-1}+\gamma \mathrm{\Delta }{t}_i{\mathbf{a}}_i$$E10

where $\beta$ and $\gamma$ are what is termed Newmark’s parameters. As formulas.

${\mathbf{u}}_i={\mathbf{u}}_{i-1}+\mathrm{\Delta }{\mathbf{u}}_i$, ${\mathbf{v}}_i={\mathbf{v}}_{i-1}+\mathrm{\Delta }{\mathbf{v}}_i$, ${\mathbf{a}}_i={\mathbf{a}}_{i-1}+\mathrm{\Delta }{\mathbf{a}}_i$ (11)

holds we can write the following formula for the vector of acceleration increments and velocity increments:

$\mathrm{\Delta }{\mathbf{a}}_i=\mathrm{\Delta }{\overline{\mathbf{a}}}_i+\frac{1}{\beta h^2}\mathrm{\Delta }{\mathbf{u}}_i$, $\mathrm{\Delta }{\mathbf{v}}_i=\mathrm{\Delta }{\overline{\mathbf{v}}}_i+\frac{\gamma }{\mathit{\beta h}}\mathrm{\Delta }{\mathbf{u}}_i$ (12)where.

$\mathrm{\Delta }{\overline{\mathbf{a}}}_i=\left(-\frac{1}{\beta \mathrm{\Delta }{t}_i}{\mathbf{v}}_{i-1}-\frac{1}{2\beta }{\mathbf{a}}_{i-1}\right)$, $\mathrm{\Delta }{\overline{\mathbf{v}}}_i=\left(1-\frac{\gamma }{2\beta }\right)\mathrm{\Delta }{t}_i{\mathbf{a}}_{i-1}-\frac{\gamma }{\beta }{\mathbf{v}}_{i-1}$ (13)

The total increments of the vector of displacement can be written as:

$$\mathrm{\Delta }{\mathbf{u}}_i=\sum _{k=1}^{n_{\mathit{\text{ITER}}}}{}^k\mathrm{\Delta }\mathrm{\Delta }{\mathbf{u}}_i$$E11

By substitution according to formulas above into the equation of motion (2) and modifying the obtained relation, we get:

$$\left(\frac{1}{\beta \mathrm{\Delta }{t_i}^2}\mathbf{M}+\frac{\gamma }{\beta \mathrm{\Delta }{t}_i}\mathbf{C}+{{}^k\mathbf{K}}_{\mathrm{T},i}\right)\cdot {}^k\mathrm{\Delta }{\mathbf{u}}_i={}^k\mathrm{\Delta }{\mathbf{f}}_i-\mathbf{M}\mathrm{\Delta }{\overline{\mathbf{a}}}_i-\mathbf{C}\mathrm{\Delta }{\overline{\mathbf{v}}}_i$$E12

where ${{}^k\mathbf{K}}_{\mathrm{T},i}$ is the tangent stiffness matrix for the kth iteration of the ith time step.

The bracket on the left-hand side of (15) represents what is termed a modified stiffness matrix, which can be denoted ${\widehat{\mathbf{K}}}_i$. The presented formula can be written in a similar form as:

$${}^k\mathrm{\Delta }\mathrm{\Delta }{\mathbf{u}}_i:\left\{\begin{array}{l}{\widehat{\mathbf{K}}}_i{}^k\mathrm{\Delta }\mathrm{\Delta }{\mathbf{u}}_i={}^k\mathrm{\Delta }\mathrm{\Delta }{\mathbf{f}}_i-\mathbf{M}\mathrm{\Delta }{\overline{\mathbf{a}}}_i-\mathbf{C}\mathrm{\Delta }{\overline{\mathbf{v}}}_i;\text{for}k=1\\ {\widehat{\mathbf{K}}}_i{}^k\mathrm{\Delta }\mathrm{\Delta }{\mathbf{u}}_i={}^k\mathrm{\Delta }\mathrm{\Delta }{\mathbf{f}}_{\mathbf{i}};\text{for}k>1\end{array}\right.$$E13

Eq. (16) enable us to calculate partial increments of displacement.

2.3 Comparison of the explicit and implicit method in seismic analysis

For this numerical study, a wall of six-floor building was used. The structure was subjected to seismic load due to the accelerogram from Umbro-Marchigiana, Italy.

This sub-chapter aims to decide what numerical method of the transient analysis of buildings exposed to earthquake is most suitable from the point of view of accuracy and processing time. The earthquake should be strong enough to cause a highly nonlinear response of the analyzed building so the modal analysis cannot be used and the usage of a numerical method for direct time integration of the equation of motion is justified.

It is widely known that the explicit method (see [7, 8, 9]) is proper for analyzing processes with very short duration, such as explosions or the collisions of vehicles, i.e., processes generally studied via transient dynamic analysis.

Because the explicit method is conditionally stable, the time step must fulfill the inequality (6). Only in this case is it possible to use the diagonal form of matrices $\mathbf{M},\mathbf{C},\mathbf{K}$ to enable conditions for maximal performance. Small time steps imposed by the relation (6) are not a substantial disadvantage in the case of earthquake analysis as the earthquake accelerograms also require very short time steps.

Implicit methods of solving a set of differential equations (see [7, 10, 11, 12]) are characterized by the fact that a system of linear equations must be solved at each time step. These methods do not require such a short time step as is needed in an explicit method. Therefore, the implicit methods are mostly used for the solution of dynamical analyses of a duration longer than several seconds.

To choose what numerical method is the most suitable for transient seismic analysis of a typical building structure (see Figure 1), a real accelerogram from Italy (see Figure 2) was used along with the Drucker-Prager material model. The accelerogram of quite strong earthquake was chosen to be strong enough to cause huge nonlinear behavior with damage and plastic yielding, where spectral analysis cannot be applied and a transient analysis using the direct integration of the equation of motion has to be used. The comparison of the explicit and implicit methods method was focused on accuracy and computational performance. A time step of 0.0001 s was used for the explicit method due to stability requirements. The same time step was also chosen for the implicit method to compare the time of processing of both methods. Another reason was to obtain a very accurate reference solution for accuracy comparisons.

Figure 1.

The analyzed concrete wall with 21 GPa of initial elasticity modulus, elements size of 1 meter, thickness of 0.2 meters [2].

Figure 2.

Accelerograms Umbro-Marchigiana, station Colfiorita-Casermette [2].

The following graphs (Figures 3–6) show the time course of the horizontal displacement of the upper right corner of the building. The graphs show very good concordance between the explicit and implicit method. It can be also seen that RFEM and ANSYS programs provide practically identical results.

Figure 3.

Explicit method—The time course of the horizontal displacement of the upper right corner for the time step 0.0001 s [2].

Figure 4.

Implicit method—The time course of the horizontal displacement of the upper right corner for the time step 0.0001 s [2].

Figure 5.

Implicit method—The time course of the horizontal displacement of the upper right corner for the time step 0.002 s [2].

Figure 6.

Implicit method—The time course of the horizontal displacement of the upper right corner for the time step 0.005 s [2].

The calculation performed by the implicit method for the same time step was approximately five times slower than the explicit method. Increasing the time step for the implicit method showed that this method is accurate enough until the time step was twenty times greater than that for the explicit method.

The implicit method provides good results for a computational time four times lower than that of the explicit method. But because of the difficulty to estimate the highest acceptable time step from the point of view of accuracy, it can be concluded that explicit and implicit methods are comparable for seismic analysis, but the implicit Newmark’s method is slightly preferred.

An interesting partial observation of this study is the finding that both basic numerical methods for the transient analysis, namely the explicit method and Newmark’s implicit method, are competitive for seismic analyses.

Results of the numerical analysis showed excellent concordance between the results of the ANSYS and RFEM programs for the implicit method. The comparison in explicit methods showed a bit worse concordance between the above-mentioned programs. The RFEM program gives more accurate results than ANSYS.

3. Plasticity and damage response of structure

The aim of the article is to compare responses of a building structure for two different nonlinear material models. A wall of a six-story building (see Figure 1) exposed to seismic load according to an accelerogram from Umbro-Marchigiana, Italy, was used to compare the two material models.

Three different material models were used: the linear elastic model, the Drucker-Prager plasticity model, and the Mazars continuous damage model [4]. Geometric nonlinearity was used for all the calculations but its influence is fairly small due to the relatively small displacement and rotation values involved.

3.1 Elasto-plastic material model

The Drucker-Prager model with isotropic hardening was executed to include different behavior of material in tension and compression which is necessary for concrete. The yield surface has the following form [3]:

$$\mathrm{\Phi }\left(\mathbf{\sigma },{\overline{\epsilon }}^p\right)=\sqrt{J_2\left(\mathbf{\sigma }\right)}+c_1\frac{1}{3}{I}_1\left(\mathbf{\sigma }\right)-c_2\mathit{coh}\left({\overline{\epsilon }}^p\right)$$E14

where $J_2$ is the second deviatoric stress invariant, $I_1$ is the first stress invariant, $\mathit{coh}$ is the cohesion that is dependent on accumulated plastic strain ${\overline{\epsilon }}^p$, and the relevant coefficients $c_1$, $c_2$ are chosen according to the required approximation to the Mohr-Coulomb criterion or can be calculated using the stress–strain diagram for uniaxial stress state and determining two conditions for the relevant yield points in tension and compression.

A standard procedure for elastic prediction and plastic correction with an implicit algorithm was used for the analysis of this elasto-plastic model. After the iterative calculation (general return-mapping update formula) of the plastic multiplicator $\mathrm{\Delta }\gamma$, the resultant stress tensor can be calculated [3]:

$${\mathbf{\sigma }}_{n+1}={\mathbf{\sigma }}_{n+1}^{\mathit{\text{trial}}}-\mathrm{\Delta }\gamma {\mathbf{D}}^e:{\mathbf{N}}_{n+1}$$E15

In this case, the updated stress, ${\mathbf{\sigma }}_{n+1}$, obtained by the implicit return mapping is the projection of the trial stress ${\mathbf{\sigma }}_{n+1}^{\mathit{\text{trial}}}$ onto the updated yield surface along the direction of the tensor ${\mathbf{D}}^e:{\mathbf{N}}_{n+1}$. Note that, since the definition of the flow vector in the smooth portion of the cone differs from that at the apex singularity, two possible explicit forms exist for the return-mapping algorithm.

3.2 Elasto-damage material model

To account for the different nonlinear performances of concrete under tension and compression, and to explicitly reproduce the dissimilar effects of the tensile and shear damage mechanisms, a decomposition of the effective stress tensor into positive and negative components with the fourth-order projection tensors ${\mathbf{P}}^t$ and ${\mathbf{P}}^c$. The Mazars isotropic damage model is modified to consider the various tensile and compressive behaviors of concrete in a more effective manner. The calculation of the resultant stress tensor is carried out with the aid of linear elastic stress estimates ${\mathbf{\sigma }}^{\mathit{\text{trial}}}$ and the damage parameters $d^t$ and $d^c$ in the following way [4]:

$$\begin{array}{l}\mathbf{\sigma }=\left(1-d^t\right){\mathbf{\sigma }}^t+\left(1-d^c\right){\mathbf{\sigma }}^c=\left(1-d^t\right){\mathbf{P}}^t:{\mathbf{\sigma }}^{\mathit{\text{trial}}}+\left(1-d^c\right){\mathbf{P}}^c:{\mathbf{\sigma }}^{\mathit{\text{trial}}}\\ =\left[\left(1-d^t\right){\mathbf{P}}^t+\left(1-d^c\right){\mathbf{P}}^c\right]:{\mathbf{\sigma }}^{\mathit{\text{trial}}}\end{array}$$E16

where ${\mathbf{\sigma }}^{\mathit{\text{trial}}}={\mathbf{D}}^{\mathbf{e}}:\mathbf{\epsilon }$.

Damage parameters are derived from the stress–strain diagram where equivalent strains calculated according to [4] are used:

$${\epsilon }_t=\frac{I_{\epsilon }}{2\left(1-2\nu \right)}+\frac{\sqrt{J_{\epsilon }}}{2\left(1+\nu \right)},$$E17

$${\epsilon }_c=\frac{I_{\epsilon }}{5\left(1-2\nu \right)}+\frac{6\sqrt{J_{\epsilon }}}{5\left(1+\nu \right)}$$E18

3.3 Comparison of the material models

The Drucker-Prager plasticity model parameters were selected so as to ensure that there is an adequate level of plasticization in the most exposed areas of the structure for the given accelerogram. The aim of the investigation is not to execute the calculations required for a real building, but merely to compare the suitability for use of various numerical methods, as well as the behavior of various material models. While the earthquake was taking place, plasticization occurred in some parts of the structure. This demonstrated itself through increased damping and a change in the distribution of tension with regard to the linear material. It can be seen from the results that at a time of 10 s (i.e., practically straight after seismicity ended) the shape of the building had undergone permanent deformation. This was in accordance with expectations. There was an overall increase in the height of the building, as well as in its width in the plasticized area. On higher floors, where plasticization did not occur, the width of the building remained the same. The increase in height and the partial expansion in terms of width might seem surprising at first sight, but these phenomena can be explained by the fact that only plasticization in tension occurred, and that the rocking of the building from one side to the other caused plasticization to take place on both sides of the structure. This increased the height of the whole building, and not only in its tilting. It also became wider for similar reasons.

The Mazars damage model expects the structure to fail due to micro cracks. The cracks are not localized but it is expected that they occur continuously (smeared cracking model). According to the Mazars model, it can be concluded from the stress–strain curve for the material that micro-cracks close after stress disappears, and that deformation also disappears in a similar way, as took place in the case of the linear elastic model. However, failure did occur. It demonstrated itself via the fact that the stress–strain curve followed a different trajectory during unloading and thus energy dissipation occurred, which resulted in a damping in oscillation. After the end of seismicity, no permanent deformation or strain remained, unlike in the case of plastic material, but there was a permanent decrease in the stiffness of the material and a change in the response of the structure to loading as a result (in this case, the loading was from continuing seismic effects). The decrease in stiffness will have an impact on dynamics in terms of a decrease in natural oscillation frequencies. The following images show the differences in the response of the structure to the same seismic loads for various material models. The corresponding results of both material models have been placed next to each other, e.g., resultant deformations and strain, with the same being shown for the time of maximum displacement (see Figure 7).

Figure 7.

The resulting shapes of the building a) Drucker-Prager plasticity model and b) the Mazars damage model.

The main aim of the study was to compare the response to seismic loading of various material models, specifically the Drucker-Prager plasticity model and Mazars damage model. At first sight, the results appeared somewhat surprising but on closer investigation, it became clear that they are correct, and that they correspond to the relevant material models. In the case of Drucker-Prager, the building remained slightly tilted after the end of the earthquake (see Figure 8, which was expected) and was also higher and wider as a whole (which was surprising at first). This somewhat remarkable result was obtained because plasticization took place only in tension and, because of the oscillation of the building from one side to the other, plasticization due to tension gradually occurred on both sides of the building. The result was an overall increase in the height of the building. For similar reasons, the building was also wider after the earthquake. In the case of the Mazars model, the geometry of the building remained the same after the end of the earthquake even though material failure occurred. However, the building was damaged. The micro-cracks did close after loading, though, and the deformation and stress from seismicity remained zero. The damage only demonstrated itself through lowered stiffness of the material in the damaged areas, which is also demonstrated in Figure 9 by the decreasing frequency. The energy loss in the form of warmth appeared as increased damping in the numerical solution.

Figure 8.

Drucker-Prager plasticity model—Time diagram of the horizontal displacement of the highest and the lowest nodes.

Figure 9.

Mazars damage model—Time diagram of the horizontal displacement of the highest and the lowest nodes.

4. Dynamic damping—Comparison of different concepts from the point of view of their physical nature and effects on civil engineering structures

There are a variety of sources of dynamical damping of vibration in civil engineering structures. Generally, damping can be caused either by the external environment or by energy dissipation due to structural deformation. The most common way to account for the damping in the motion equation is Rayleigh damping (see [13, 14, 15]). The Rayleigh damping is proportional to a linear combination of the stiffness and mass matrix.

This subchapter deals with the justification for the use of Rayleigh damping and discusses another solution because there is no direct physical interpretation of the mass damping parameter and the stiffness damping parameter. The damping given by internal resistance of structure occurs in the case of using inelastic materials when the loading and unloading parts of the strain–stress diagram differ, and therefore dissipation occurs in loading cycles, as described in [1, 2, 3, 5, 6, 10, 11, 12, 16, 17, 18, 19, 20]. Another source of damping is friction in structural connections.

4.1 Finite element analysis of the equation of motion

The virtual work of external forces has to be equal to the virtual work of internal forces. For one element with the volume ${\mathrm{\Omega }}_e$ and the surface boundary ${\mathrm{\Gamma }}_e$ we can write this equilibrium of virtual work as follows:

$$\underset{{\mathrm{\Omega }}_e}{\int }\delta {\mathbf{u}}^T\mathbf{b}d{\mathrm{\Omega }}_e+\underset{{\mathrm{\Gamma }}_e}{\int }\delta {\mathbf{u}}^T\mathbf{t}d{\mathrm{\Gamma }}_e+\sum _{i=1}^n\delta {\mathbf{u}}_i^T{\mathbf{f}}_i=\underset{{\mathrm{\Omega }}_e}{\int }\left(\delta {\mathbf{u}}^T\rho \ddot{\mathbf{u}}+\delta {\mathbf{u}}^{T}c\dot{\mathbf{u}}+\delta {\mathbf{\epsilon }}^T\mathbf{\sigma }\right)d{\mathrm{\Omega }}_e$$E19

where $\delta \mathbf{u}$ and $\delta \mathbf{\epsilon }$ represent the virtual displacement and pertinent strain, respectively, $\rho$ is the material density, $c$ is the viscous damping parameter, $\mathbf{b}$ and $\mathbf{t}$ are the volume and surface forces respectively, ${\mathbf{f}}_i$ and $\delta {\mathbf{u}}_i$ represent the concentrated forces and pertinent generalized displacement, respectively. The damping parameter $c$ does not take into account material (internal) or external resistance. The material viscosity is not taken into account because it has an impact only in the case of a nonzero rate of deformation in the given mass point (strain of the body). The influence of the external environment is not taken into account because the stress vector $c\dot{\mathbf{u}}$, which acts against the motion $\mathbf{u}$, arises even in the internal points of the bodies.

When discretizing for the finite element method we obtain the following relations:

$$\mathbf{u}=\mathbf{Nd}\dot{\mathbf{u}}=\mathbf{N}\dot{\mathbf{d}}\ddot{\mathbf{u}}=\mathbf{N}\ddot{\mathbf{d}}\mathbf{\epsilon }=\mathbf{Bd}$$E20

Combining eqs. (22) and (23) we obtain:

$$\delta {\mathbf{d}}^T\left(\underset{{\mathrm{\Omega }}_e}{\int }\rho {\mathbf{N}}^T\mathbf{N}d{\mathrm{\Omega }}_e\ddot{\mathbf{d}}+\underset{{\mathrm{\Omega }}_e}{\int }c{\mathbf{N}}^T\mathbf{N}d{\mathrm{\Omega }}_e\dot{\mathbf{d}}+\underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T\mathbf{\sigma }d{\mathrm{\Omega }}_e-\underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{N}}^T\mathbf{b}d{\mathrm{\Omega }}_e-\underset{{\mathrm{\Gamma }}_e}{\int }{\mathbf{N}}^T\mathbf{t}d{\mathrm{\Gamma }}_e-\sum _{i=1}^n{\mathbf{f}}_i\right)=0$$E21

where it is assumed that the concentrated forces ${\mathbf{f}}_i$ act in nodes. Let us denote the first two integrals in the equation as the consistent mass matrix and the damping matrix.

$${\mathbf{M}}_e=\underset{{\mathrm{\Omega }}_e}{\int }\rho {\mathbf{N}}^T\mathbf{N}d{\mathrm{\Omega }}_e$$E22

$${\mathbf{C}}_e=\underset{{\mathrm{\Omega }}_e}{\int }c{\mathbf{N}}^T\mathbf{N}d{\mathrm{\Omega }}_e$$E23

The word consistent means that the matrix follows directly from the discretization of a finite element with corresponding shape functions $\mathbf{N}$. Let us define the vectors of the internal and external nodal forces:

$${\mathbf{f}}_e^{\mathrm{int}}=\underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T\mathbf{\sigma }d{\mathrm{\Omega }}_e$$E24

$${\mathbf{f}}_e^{\mathit{ext}}=\underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{N}}^T\mathbf{b}d{\mathrm{\Omega }}_e+\underset{{\mathrm{\Gamma }}_e}{\int }{\mathbf{N}}^T\mathbf{t}d{\mathrm{\Gamma }}_e+\sum _{i=1}^n{\mathbf{f}}_i$$E25

Substituting from eqs. (25), (26), (27), and (28) into eq. (24), and taking into account the fact that variation $\delta \mathbf{d}$ can be arbitrary, and hence the second form of the product must be equal to zero, we obtain:

$${\mathbf{M}}_e\ddot{\mathbf{d}}+{\mathbf{C}}_e\dot{\mathbf{d}}+{\mathbf{f}}_e^{\mathrm{int}}={\mathbf{f}}_e^{\mathit{ext}}$$E26

For a linear elastic material without viscosity, we can write the following relation for the internal nodal forces:

$${\mathbf{f}}_e^{\mathrm{int}}={\mathbf{K}}_e\mathbf{d}$$E27

where

$${\mathbf{K}}_e=\underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T\mathbf{DB}d{\mathrm{\Omega }}_e$$E28

is the stiffness matrix of the element, with $\mathbf{D}$ being the constitutive matrix of the material. Then eq. (29) can be rewritten as:

$${\mathbf{M}}_e\ddot{\mathbf{d}}+{\mathbf{C}}_e\dot{\mathbf{d}}+{\mathbf{K}}_e\mathbf{d}={\mathbf{f}}_e^{\mathit{ext}}$$E29

Eqs. (30) and (32) express in discretized form Newton’s second law, or, more generally, the law of conservation of momentum. When writing these equations in the global form, i.e., for all degrees of freedom of the analyzed structure, we obtain the following equation:

$$\mathbf{M}\ddot{\mathbf{d}}+\mathbf{C}\dot{\mathbf{d}}+{\mathbf{f}}^{\mathrm{int}}={\mathbf{f}}^{\mathit{ext}}$$E30

$$\mathbf{M}\ddot{\mathbf{d}}+\mathbf{C}\dot{\mathbf{d}}+\mathbf{Kd}={\mathbf{f}}^{\mathit{ext}}$$E31

4.2 Rayleigh damping

For Rayleigh damping, the damping matrix ${\mathbf{C}}_e^R$ is defined as the linear combination of the consistent mass matrix ${\mathbf{M}}_e$ and the stiffness matrix ${\mathbf{K}}_e$.

$${\mathbf{C}}_e^R=\alpha {\mathbf{M}}_e+\beta {\mathbf{K}}_e$$E32

When substituting to eq. (35) ${\mathbf{M}}_e$ from eq. (25) and ${\mathbf{K}}_e$ from (31), we obtain the relation for the damping matrix in the following form:

$${\mathbf{C}}_e^R=\alpha \underset{{\mathrm{\Omega }}_e}{\int }\rho {\mathbf{N}}^T\mathbf{N}d{\mathrm{\Omega }}_e+\beta \underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T\mathbf{DB}d{\mathrm{\Omega }}_e$$E33

The damping matrix has two parts and the first one is identical with ${\mathbf{C}}_e$ as defined in (26), where $c=\mathit{\alpha \rho }$ and its physical nature is unclear. The second part of the expression does not correspond with relation (26) but is proportional to the stiffness matrix. The damping for $\alpha =0$ and $\beta >0$ is thus proportional to the rate of deformation of the body. The internal nodal forces arise only when the body deforms over time and no internal forces arise when an element moves as a rigid body.

4.2.1 Damping caused by material viscosity

As we mentioned before the internal damping caused by material arises only when a body is strained and is proportional to the strain rate. We have chosen the three most well-known viscoelastic models for this article for the damping caused by material viscosity and they are described below.

4.2.2 The Maxwell material model

For the Maxwell (see Figure 10) model the following relations are valid:

Figure 10.

Scheme of the Maxwell model.

$$\sigma \left(t\right)={\sigma }^e\left(t\right)={\sigma }^v\left(t\right),\epsilon \left(t\right)={\epsilon }^e\left(t\right)+{\epsilon }^v\left(t\right),{\epsilon }^e\left(t\right)=\frac{\sigma \left(t\right)}{E},{\dot{\epsilon }}^v\left(t\right)=\frac{\sigma \left(t\right)}{\eta }$$E34

As a result, we obtain the following linear nonhomogenous ordinary differential equation which describes the constitutive relation between stress and strain:

$$\sigma \left(t\right)+\frac{\eta }{E}\dot{\sigma }\left(t\right)=\eta \dot{\epsilon }\left(t\right)$$E35

4.2.3 The kelvin-Voigt material model

The Kelvin-Voigt model (see Figure 11) is based on the following relations:

Figure 11.

Scheme of the kelvin-Voigt model.

$$\epsilon \left(t\right)={\epsilon }^e\left(t\right)={\epsilon }^v\left(t\right),\sigma \left(t\right)={\sigma }^e\left(t\right)+{\sigma }^v\left(t\right),{\epsilon }^e\left(t\right)=\frac{\sigma \left(t\right)}{E},{\dot{\epsilon }}^v\left(t\right)=\frac{\sigma \left(t\right)}{\eta }$$E36

Using these relations, we obtain as a result the following linear non-homogenous ordinary differential equation which describes the constitutive relation between stress and strain:

$$\sigma \left(t\right)=\mathit{E\epsilon }\left(t\right)+\eta \dot{\epsilon }\left(t\right)$$E37

4.2.4 The standard linear solid (SLS) material model

For the SLS model (see Figure 12) the following relations are valid:

Figure 12.

Scheme of the standard linear solid model.

$$\epsilon \left(t\right)={\epsilon }_1\left(t\right)={\epsilon }_2\left(t\right),{\epsilon }_1\left(t\right)={\epsilon }_1^e\left(t\right)=\frac{\sigma \left(t\right)}{E_1},{\epsilon }_2\left(t\right)={\epsilon }_2^e\left(t\right)+{\epsilon }_2^v\left(t\right),\sigma \left(t\right)={\sigma }_1\left(t\right)+{\sigma }_2\left(t\right),$$

$${\sigma }_1\left(t\right)=E_1\epsilon \left(t\right)and{\sigma }_2\left(t\right)=E_2\left(\epsilon \left(t\right)-{\epsilon }_2^v\left(t\right)\right)$$E38

Then, for the resulting stress the following relations can be written:

$$\sigma \left(t\right)=E_1\epsilon \left(t\right)+E_2\left(\epsilon \left(t\right)-{\epsilon }_2^v\left(t\right)\right)$$E39

After several modifications, we obtain the final relation for expressing strain rate as a function of stress rate, stress, and actual strain (it is a differential equation describing the constitutive relation between stress and strain).

$$\sigma \left(t\right)+\frac{\eta }{E_2}\dot{\sigma }\left(t\right)=\frac{\eta \left(E_1+E_2\right)}{E_2}\dot{\epsilon }\left(t\right)+E_1\epsilon \left(t\right)$$E40

4.2.5 Damping caused by the kelvin-Voigt model

Let us now demonstrate how the damping caused by the Kelvin-Voigt model will manifest itself in the equation of motion. Let us substitute for $\mathbf{\sigma }\left(t\right)$ from (41) into eq. (27). For the vector of internal nodal forces, we obtain:

$${\mathbf{f}}_e^{\mathrm{int}}=\underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T\left({\mathbf{\sigma }}^e+\eta \dot{\mathbf{\epsilon }}\right)d{\mathrm{\Omega }}_e=\underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T\left({\mathbf{\sigma }}^e+\eta \mathbf{B}\dot{\mathbf{d}}\right)=\underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T{\mathbf{\sigma }}^e+\eta \underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T\mathbf{B}\dot{\mathbf{d}}$$E41

Taking into consideration the fact that the damping will be caused by the material, we will not expect damping to take place with the help of the matrix $\mathbf{C}$ here. The equation of motion shall then read:

$${\mathbf{M}}_e\ddot{\mathbf{d}}+\underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T{\mathbf{\sigma }}^e+\eta \underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T\mathbf{B}\dot{\mathbf{d}}={\mathbf{f}}_e^{\mathit{ext}}$$E42

The equation can then be rewritten for a linear elastic material as:

$${\mathbf{f}}_e^{\mathrm{int}}={\mathbf{K}}_e\mathbf{d}+\eta \underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T\mathbf{B}\dot{\mathbf{d}}d{\mathrm{\Omega }}_e$$E43

$${\mathbf{M}}_e\ddot{\mathbf{d}}+{\mathbf{f}}_e^{\mathrm{int}}={\mathbf{f}}_e^{\mathit{ext}}$$E44

and after the substitution for ${\mathbf{f}}_e^{\mathrm{int}}$, we can rewrite the equation of motion in the form

$${\mathbf{M}}_e\ddot{\mathbf{d}}+{\mathbf{K}}_e\mathbf{d}+\eta \underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T\mathbf{B}\dot{\mathbf{d}}d{\mathrm{\Omega }}_e={\mathbf{f}}_e^{\mathit{ext}}$$E45

4.3 Comparison of Rayleigh damping with the damping of the kelvin-Voigt material

The Rayleigh damping depends only on the velocities of mass points in space, and that damping also arises in the case of the movement of a rigid body due to coefficient $\alpha$. To compare the Rayleigh damping and damping caused by the Kelvin-Voigt material model let us assume $\alpha =0$. The equation of motion then will read as:

$${\mathbf{M}}_e\ddot{\mathbf{d}}+\beta \underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T\mathbf{DB}d{\mathrm{\Omega }}_e\dot{\mathbf{d}}+\underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T\mathbf{\sigma }d{\mathrm{\Omega }}_e={\mathbf{f}}_e^{\mathit{ext}}$$E46

If we compare the equation with the Kelvin-Voigt material model (45) with this eq. (49), one can see that the difference between them lies only in the term with the first derivative of the deformation parameters by time, i.e., the velocities $\dot{\mathbf{d}}$. If we simplify (reduce) these equations to 1D tasks, substituting Young’s modulus $E$ for the constitutive matrix $\mathbf{D}$ and comparing the terms with $\dot{\mathbf{d}}$ from the equations we obtain the relation

$$\mathit{\beta E}\underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T\mathbf{B}d{\mathrm{\Omega }}_e\dot{\mathbf{d}}=\eta \underset{{\mathrm{\Omega }}_e}{\int }{\mathbf{B}}^T\mathbf{B}d{\mathrm{\Omega }}_e\dot{\mathbf{d}}$$E47

and subsequently simple relation between the parameters $\beta$ and $\eta$ of both models.

$$\boxed{\mathit{\beta E}=\eta }$$E48

It is clear that the Rayleigh damping and the damping caused by Kelvin-Voight material model are identical for 1D problems.

4.4 Damping caused by the plasticizing, or damaging of material

Figure 13 shows a loading and unloading diagram for a) elasto-plastic b) elasto-damage material for the 1D stress and strain state.

Figure 13.

Loading and unloading diagram for elasto-plastic (a) and elasto-damage (b) materials, with division of energy into elastic and dissipative.

$$\text{a)}{w}^p\left(\overline{t}\right)=\underset{o}{\overset{\overline{t}}{\int }}\mathbf{\sigma }\left(t\right):{\dot{\mathbf{\epsilon }}}^p\left(t\right)\mathit{dt}=\underset{o}{\overset{\overline{t}}{\int }}\mathbf{\sigma }\left(t\right):\dot{\mathbf{\epsilon }}\left(t\right)\mathit{dt}-\frac{1}{2}{\mathbf{C}}^{-1}:\overline{\mathbf{\sigma }}\left(\overline{t}\right):\overline{\mathbf{\sigma }}\left(\overline{t}\right)$$E49

$$\text{b)}{w}^p\left(\overline{t}\right)=\underset{o}{\overset{\overline{t}}{\int }}\mathbf{\sigma }\left(t\right):{\dot{\mathbf{\epsilon }}}^p\left(t\right)\mathit{dt}=\underset{o}{\overset{\overline{t}}{\int }}\mathbf{\sigma }\left(t\right):\dot{\mathbf{\epsilon }}\left(t\right)\mathit{dt}-\frac{1}{2}\overline{\mathbf{\epsilon }}\left(\overline{t}\right):\overline{\mathbf{\sigma }}\left(\overline{t}\right)$$E50

This dissipation energy is lost from the mechanical system in the form of heat and this process manifests itself in dynamics as damping. The sum of this dissipation energy $w^p$ plus the elastic energy $w^e$ is equal to the total energy used for the deformation work.

In the case of plasticizing the damping occurs at the time when a plastic deformation, arises, and later vibration is no longer damped unless the maximum strain achieved so far is overcome once again. Plasticizing does not affect the natural frequency of the structure.

In the case of damage, the damping occurs at the time when a material damage arises, and later vibration is no longer damped unless the maximum strain achieved so far is overcome once again. Damage causes material softening, so the Eigen frequencies will be decreased.

4.5 Damping caused by friction

The Coulomb friction (see the diagram in Figure 14) is also a significant part of internal damping of structures. It occurs in connections of structural elements and most often comes about in screw or rivet connections in steel structures. Mechanical work and pertinent dissipation arise in the case of any relative motion in the connections. In contrast to plastic or damage behavior, mechanical work is performed with each relative forward and backward motion, and the pertinent damping is then permanent during vibration. The friction force is usually proportional to the pressure force in connections (Coulomb friction).

Figure 14.

Friction diagram [1].

5. Conclusions

The conclusion is that for accurate nonlinear seismic analysis we cannot use spectral methods but we have to use the methods of direct integration of the motion equation. We can use both methods (explicit and implicit) with a slight advantage of the explicit method with respect to time consumption. Another important conclusion is that the mathematical model of the concrete has a significant influence on the nonlinear transient response of structure (earthquake analysis) and reality is somewhere between analyzed elasto-plastic and elasto-damage behavior and thus we recommend using combination of plasticity and damage models simultaneously. The analysis of different sources of damping showed that it is needed to use not only external damping but also (and mainly) internal damping using material inelasticity models (viscosity, plasticity/damage, friction) instead of simple Rayleigh damping which is insufficient for the huge nonlinear time-dependent response of structure on earthquake loading.

Acknowledgments

The work presented in this Chapter has been supported by the project of specific university research at Brno University of Technology No. FAST-S-22-7867.