Time history analysis formulation in SCAD FEA software

Sergiy Fialko1 , Viktor Karpilovskyi2

1Tadeusz Kościuszko Cracow University of Technology, Cracow, Poland

2IT Company SCAD Soft, Kiev, Ukraine

1Corresponding author

Journal of Measurements in Engineering, Vol. 6, Issue 4, 2018, p. 173-180. https://doi.org/10.21595/jme.2018.20408
Received 16 October 2018; received in revised form 23 November 2018; accepted 30 November 2018; published 31 December 2018

Copyright © 2018 Sergiy Fialko, et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License

This paper provides the formulation of the problem of forced vibrations of structures, which has the following peculiarities: equations of motion are formulated in absolute coordinates allowing to take into account the asynchronous excitations of the supports, smoothing of the time functions for forced displacements with the help of Hermite polynomials, and also the possibility of considering damping that does not obey the Rayleigh hypothesis.

Keywords: finite element method, structural dynamics, seismic analysis, asynchronous excitations, local and material damping.

1. Introduction

In the classical formulation of the seismic problem, the motion of the finite element design model is considered in the coordinate system associated with the mobile rigid platform supporting the entire structure (Fig. 1).

Fig. 1. The classical formulation of a seismic problem

Fig. 2. The formulation of a seismic problem in absolute coordinates

The equations of motion are written as follows:

M u ¨ r + C u ˙ r + K u r = - M u ¨ e t ,

here M, C, K are mass, dissipation and stiffness matrices respectively, u is the displacement vector, the subscript e denotes the components of the bulk motion – movement of the platform points with respect to a fixed coordinate system Oxyz. Components of the relative motion – movement of the points of the structure with respect to the mobile platform, denoted by the subscript r, are the unknowns. The input action is the given accelerogram u¨e(t).

In the case of the asynchronous excitation of the supports the equations of motion are written as follows [1]:

M u ¨ d + C u ˙ d + K u d = -   M K - 1 K 12 + M 12 u ¨ e t ,
u t = u d t + u s t ,

where ut are the displacements of nodes of the design model in the absolute coordinate system fixed to the stationary ground, udt are the displacements with respect to the quasistatic state ust, determined by solving the problem:

K u s = - K 12 u ¨ e t .

In Eqs. (2) and (4) K12 and M12 are stiffness and mass submatrices corresponding to the degrees of freedom for which the accelerations u¨e(t) are given. As a rule, the problems Eqs. (1) and (2) are solved by the method of expansion in natural modes determined for a structure with fixed supports (ue= 0). In order to separate the equations of motion Eqs. (1), (2) into unrelated equations in normal coordinates, the dissipation matrix C has to be represented as a linear combination of the mass and stiffness matrices, i.e. the Rayleigh hypothesis has to be fulfilled. However, in the case of multi-component damping, when different structural elements are made of materials with different damping, and also in the case of local dampers, the application of this approach causes significant difficulties.

Moreover, the approach Eq. (2) in the case of an asynchronous excitation of supports requires an additional solution of the quasistatic problem Eq. (4) with subsequent superposition Eq. (3) of solutions of the dynamic and static problems, which substantially complicates the solution procedure.

The equations of motion are written in absolute coordinates in SCAD [2], (Fig. 2), which provides the natural solution of the problem of forced vibrations at the asynchronous excitation of supports both in the case of the Rayleigh damping and the multi-component damping, and in the case of the local dampers as well. The proposed formulation of the problem is presented in this paper.

2. Problem formulation in absolute coordinates

The equation of motion (Fig. 2) has the form:

0 0 0 M u ¨ e u ¨ + C 11 C 12 C 21 C u ˙ e u ˙ + K 11 K 12 K 21 K u e u = 0 p t ,

where the subscript e corresponds to the specified displacements, velocities and accelerations of the bulk motion, and pt is the vector of the external forces of the finite element design model. In the case of the seismic action pt= 0. It follows from Eq. (5) that:

M u ¨ + C u ˙ + K u = - K 12 u e t - C 12 u ˙ e t + p t   .

The left-hand side Eq. (6) contains displacements, velocities and accelerations that have to be determined, and the right-hand side expression contains known displacements and velocities. All components of the displacement, velocity and acceleration vectors are given in the absolute coordinate system Oxyz fixed to the stationary ground. Dissipation matrix blocks C and C12 are assembled in exactly the same way as the stiffness matrix blocks K and K12, using a standard assembly procedure covering the finite elements of the design model [3]. It is stated in [1] that the term C12u˙et in Eq. (6) can be neglected. However, we will keep this term for the sake of precision. In the proposed approach the input signal is not an accelerogram, as is in the case of the classical approach described by Eqs. (1) and (2), but a vibrogram uet and a velocigram u˙et. Therefore, the given accelerogram u¨et must be integrated twice, using the Accelerogram Editor from SCAD Office [2].

If Eq. (6) is complemented with the initial conditions (usually accepted: u0= u˙0= 0), we obtain the Cauchy problem. The direct implicit α – HHT method [4] is usually used for its numerical integration.

3. Smoothing of forced displacements

Let’s consider the problem of vibrations of a linear oscillator excited by a given displacement of the support (Fig. 3).

Fig. 3. The linear oscillator with a rigid link

Fig. 4. The given acceleration and displacement in the node 1

A peculiarity of this problem is modeling a rigid link connecting nodes 1 and 2 using the method of penalty functions [5]. This approach is widely used in SCAD and has a number of advantages over the traditional method of accounting for rigid links, based on kinematic and static relationships of structural mechanics. The physical essence of the method of penalty functions is equivalent to setting springs, the stiffness of which is several orders of magnitude greater than the stiffness of the finite elements adjacent to the nodes connected by a rigid link. In the exact solution of this problem, the horizontal displacements, velocities, and accelerations of the node 2 must be exactly equal to the respective displacements, velocities, and accelerations of the node 1.

The Newmark method is used for the numerical solution (the α parameter in the α – HHT method is taken as zero). Displacements of the nodes 1 and 2 turn out to be practically equal, the velocities differ insignificantly, and the horizontal acceleration of the node 2 is completely distorted by the high-frequency component (Fig. 5), which appeared due to the high stiffness caused by the application of the penalty function method in the modeling of a rigid link. It should be noted that the displacements, velocities and accelerations of the node 3, as well as the forces and bending moments in the bar are obtained correctly.

Fig. 5. Acceleration in the node 2. High-frequency component distorts the solution

The reason for the distortion of the acceleration in the node 2 was the application of the linear interpolation to approximate the time function of the given displacements of the node 1. The velocity-time relationship is represented by a piecewise constant function with discontinuities at the interpolation nodes (Fig. 6). Accelerations between the interpolation nodes are zero, and those at the nodes are represented by Dirac delta functions.

Fig. 6. Approximation of the displacements, velocities, and accelerations when using the linear interpolation for the given displacements

While the linear interpolation was used to approximate the time function for the external forces pt, the inertial properties of the design model were smoothing the jumps shown in the Fig. 6. In the classical approach to seismic problems Eqs. (1) and (2), the input action is a given acceleration, so applying the linear interpolation to the time function for accelerations also does not lead to the excitation of high-frequency components.

In order to eliminate parasitic high-frequency oscillations arising in the proposed approach Eq. (6) when using forced displacements and velocities as input, we apply Hermite polynomials for the interpolation of forced displacements between the nodes:

u e i ξ = N 1 ξ     N 2 ξ     N 3 ξ     N 4 ξ d i - 1     v i - 1     d i     v i T ,

where ξ=ξi=t-ti-1/Δti,N1ξ=2ξ3-3ξ2+1,N2ξ=-Δtiξ3-2ξ2+ξ,N3ξ=-2ξ3+3ξ2, N4ξ=-Δtiξ3-ξ2, di-1, vi-1, di, vi are the displacements and velocities in the interpolation nodes i-1 and i (Fig. 7). The superscript i over the displacement u denotes the interval number, which is circled in the Fig. 7. Therefore, here ξ= ξi. It follows from Eq. (7) that:

u ˙ e i ξ =   1 t i N ˙ 1 ξ   N ˙ 2 ξ   N ˙ 3 ξ   N ˙ 4 ξ   d i - 1   v i - 1   d i   v i T ,
u ¨ e i ξ =   1 t i 2 N ¨ 1 ξ   N ¨ 2 ξ   N ¨ 3 ξ   N ¨ 4 ξ   d i - 1   v i - 1   d i   v i T .

This approach provides the absence of discontinuities of the displacements and velocities in the interpolation nodes:

u e i ξ i = 1 =   u e i + 1 ξ i + 1 = 0 ,       u ˙ e i ξ i = 1 =   u ˙ e i + 1 ξ i + 1 = 0 .

The values of the displacements in the interpolation nodes are given, and the values of the velocities are not. Our task is to select such nodal values of the velocities so that there are no discontinuities in the accelerations in the nodes:

u ¨ e i ξ i = 1 =   u ¨ e i + 1 ξ i + 1 = 0 .

Applying Eq. (11) to each internal interpolation node and supplementing it with the boundary conditions v0= vN= 0, where N is the number of the last point where the displacements are specified, we obtain a system of linear algebraic equations with a tridiagonal symmetric matrix, solving which, we determine the values of the velocities vi, i= 1, 2, …, N-1 at the interpolation nodes.

Fig. 7. Interpolation nodes and local coordinates

As a result of applying the proposed approach based on the interpolation of forced displacements by Hermite polynomials, the horizontal acceleration of the node 2 turned out to be exactly the same as the horizontal acceleration of the node 1 (Fig. 4).

4. Asynchronous excitation of supports

If a structure has a considerable length, it becomes necessary to take into account the fact that the excitations in the soil reach some supports before the others ([7], [8] and so on). Fig. 8 shows the general view of the New Safe Confinement over the fourth power unit of the Chernobyl nuclear power plant, Fig. 9 shows the design model created in SCAD, and Fig. 10 shows an accelerogram of the input action [6].

Fig. 8. The new confinement “Arka” before the sliding over the old confinement, covering the damaged structures of the fourth power unit of the Chernobyl nuclear power plant

According to the geophysical surveys [6], the velocity of the seismic wave propagation vs is 200-300 m/s, and the predominant period of ground vibration is Tprev= 0.4 s. Thus, the seismic wave length lies in the range λs=vsTprev= 80-120 m. When the seismic wave propagates along the Oy axis, the delay time of the wave arrival at the distant support is t0=L/vs= 0.86-1.29 s, i.e. t0= (2.14-3.2)∙Tprev.

The Rayleigh damping model is used. According to recommendations [6], the factor for the mass matrix is assumed to be 0.16, and the factor for the stiffness matrix is 0.01.

Fig. 9. The design model. The red arrow depicts the direction of a seismic input motion

Fig. 10. Accelerogram of the input motion

Fig. 11 shows the longitudinal force-time relationship in the element #1687 (Fig. 9), of the bottom chord of the middle arch truss. This element turned out to be one of the most loaded structural elements under the seismic action.

Fig. 11. The longitudinal force in the element #1687, t0= 0 – synchronous excitation, t0= 1.2 s – asynchronous excitation

The delay in achieving a seismic excitation of a distant support leads to a significant change in the longitudinal force-time relationship. Fig. 11 shows a time interval, where the longitudinal force reaches its largest absolute values. Although the maximum absolute values at t0= 0 and t0= 1.2 s differ insignificantly in this example, the character of these relationships are completely different, and as a result the maximum force values are reached at different moments in time.

5. Multi-component damping

Fig. 12 shows a design model of the spatial frame, the lower two floors of which are made of reinforced concrete (γ=δ/π= 0.1, γ is the coefficient of internal inelastic resistance, δ is the logarithmic decrement of oscillations), and the upper two floors are made of steel (γ= 0.03). Soil compliance in the direction of the translational displacements is modeled by the springs with stiffness kx=ky=kz= 10 MN/m, installed in the level of the supports, and the soil dissipation is modeled by local dampers with viscosity βx=βy=βz= 10 MN∙s/m. The oscillations are excited by the synchronous movement of the supports in the direction of the Oy axis, and the time function is taken as ft=1-cos(Ωt), where the frequency of the forced oscillations is Ω=ω2, and ω2 is the second natural frequency of oscillations, at which the structure moves mainly in the direction of the ground excitation.

Since the seismic action takes the structure through all the resonance modes, we set such a testing load so as to create the resonance conditions. The only factor limiting the resonance amplitude of the oscillations of the linear design model is damping. Therefore, it is very important to take into account the damping correctly under seismic actions. Curve 1 (Fig. 13) describes the horizontal displacements of the node 5 in the direction of load action when only local dampers modeling the soil dissipation are taken into account, and curve 2 describes them in the case when both the soil dissipation and damping in structural elements are taken into account.

It should be noted that a dissipation model not corresponding to the Rayleigh hypothesis was used when solving this problem. The dissipation matrix for column and beam finite elements is determined as follows:

C e = γ e K e ,

where Ce, Ke are dissipation and stiffness matrices of the finite element e, γe is the coefficient of internal inelastic resistance of the material.

Fig. 12. Design model for a reinforced concrete and steel spatial frame

Fig. 13. Horizontal displacement of the node 5. 1 – only local dampers, 2 – local dampers and material damping

6. Conclusions

The presented formulation of the problem of forced structural vibrations enables to solve problems for large span structures that require the consideration of the asynchronous excitation of supports and to consider both the Rayleigh damping model and the material dumping, that allows one on the creation of the design models with multi-component damping. It is shown that with the forced excitation of supports, linear interpolation of displacements in time can lead to the appearance of parasitic high-frequency oscillations. A method for constructing a local cubic approximation based on the use of shape functions in the form of Hermite polynomials, which does not lead to the above-mentioned high-frequency oscillations, is also proposed.


  1. Clough R. W., Penzien J. Dynamics of Structures. Computers and Structures Inc., Berkeley, CA, USA, 2003. [CrossRef]
  2. Karpilovskyi V. S., Kryksunov E. Z., Maliarenko A. A., Perelmuter A. V., Perelmuter M. A., Fialko S. Y. SCAD Office. Version 21, System SCAD++, SCAD Soft, 2018, https://scadsoft.com/download/SCAD1033.pdf. [CrossRef]
  3. Zienkiewicz O. C., Taylor R. L. The Finite Element Method. Fifth Edition, Vol. 2, Solid Mechanics, Butterworth – Heinemann, 2000. [CrossRef]
  4. Miranda I., Ferencz R. M., Hughes T. J. R. An improved implicit-explicit time integration method for structural dynamics. Earthquake Engineering and Structural Dynamics, Vol. 18, 1989, p. 643-653. [Publisher]
  5. Fialko S. Y. Application of rigid links in structural design models. International Journal for Computational Civil and Structural Engineering, Vol. 13, Issue 3, 2017, p. 119-137. [Publisher]
  6. Seismic Micro Districting of the Chernobyl NPP Industrial Site and Seismic Monitoring of the Shelter Facility. Report S. I. Subbotin Institute of Geophysics. NAS of Ukraine, Kiev, 1995-1996, (in Russian). [CrossRef]
  7. Dulińska J., Zięba A. Dynamic response of pipelines to non-uniform kinematic excitation. Technical Transactions, Civil Engineering, 2005, p. 3-21. [CrossRef]
  8. Zembaty Z. Vibrations of bridge structure under kinematic wave excitations. Journal of Structural Engineering, Vol. 123, Issue 4, 1997, p. 479-488. [Publisher]