Time history analysis formulation in SCAD FEA software

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. 1The classical formulation of a seismic problem

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

The equations of motion are written as follows:

here $\mathbf{M}$, $\mathbf{C}$, $\mathbf{K}$ are mass, dissipation and stiffness matrices respectively, $\mathbf{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 ${\ddot{\mathbf{u}}}_e\left(t\right)$.

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

where $\mathbf{u}\left(t\right)$ are the displacements of nodes of the design model in the absolute coordinate system fixed to the stationary ground, ${\mathbf{u}}_d\left(t\right)$ are the displacements with respect to the quasistatic state ${\mathbf{u}}_s\left(t\right)$, determined by solving the problem:

In Eqs. (2) and (4) ${\mathbf{K}}_{12}$ and ${\mathbf{M}}_{12}$ are stiffness and mass submatrices corresponding to the degrees of freedom for which the accelerations ${\ddot{\mathbf{u}}}_e\left(t\right)$ 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 (${\mathbf{u}}_e=$ 0). In order to separate the equations of motion Eqs. (1), (2) into unrelated equations in normal coordinates, the dissipation matrix $\mathbf{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:

where the subscript $e$ corresponds to the specified displacements, velocities and accelerations of the bulk motion, and $\mathbf{p}\left(t\right)$ is the vector of the external forces of the finite element design model. In the case of the seismic action $\mathbf{p}\left(t\right)=$ 0. It follows from Eq. (5) that:

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 $\mathbf{C}$ and ${\mathbf{C}}_{12}$ are assembled in exactly the same way as the stiffness matrix blocks $\mathbf{K}$ and ${\mathbf{K}}_{12}$, using a standard assembly procedure covering the finite elements of the design model [3]. It is stated in [1] that the term ${\mathbf{C}}_{12}{\dot{\mathbf{u}}}_e\left(t\right)$ 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 ${\mathbf{u}}_e\left(t\right)$ and a velocigram ${\dot{\mathbf{u}}}_e\left(t\right)$. Therefore, the given accelerogram ${\ddot{\mathbf{u}}}_e\left(t\right)$ must be integrated twice, using the Accelerogram Editor from SCAD Office [2].

If Eq. (6) is complemented with the initial conditions (usually accepted: $\mathbf{u}\left(0\right)=\dot{\mathbf{u}}\left(0\right)=$ 0), we obtain the Cauchy problem. The direct implicit $\alpha$ – 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. 3The linear oscillator with a rigid link

Fig. 4The 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 $\alpha$ parameter in the $\alpha$ – 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. 5Acceleration 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. 6Approximation 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 $\mathbf{p}\left(t\right)$, 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:

where $\xi ={\xi }_i=\left(t-t_{i-1}\right)/\mathrm{\Delta }{t}_i\text{,}$$N_1\left(\xi \right)=2{\xi }^3-3{\xi }^2+1\text{,}$$N_2\left(\xi \right)=-\mathrm{\Delta }{t}_i\left({\xi }^3-2{\xi }^2+\xi \right)\text{,}$$N_3\left(\xi \right)=-2{\xi }^3+3{\xi }^2$, $N_4\left(\xi \right)=-\mathrm{\Delta }{t}_i\left({\xi }^3-{\xi }^2\right)$, $d_{i-1}$, $v_{i-1}$, $d_i$, $v_i$ 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 $\xi ={\xi }_i$. It follows from Eq. (7) that:

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

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:

Applying Eq. (11) to each internal interpolation node and supplementing it with the boundary conditions $v_0=v_N=$ 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 $v_i$, $i=$ 1, 2, …, $N-1$ at the interpolation nodes.

Fig. 7Interpolation 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. 8The 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 $v_s$ is 200-300 m/s, and the predominant period of ground vibration is $T_{prev}=$ 0.4 s. Thus, the seismic wave length lies in the range ${\lambda }_s=v_s\bullet T_{prev}=$ 80-120 m. When the seismic wave propagates along the $Oy$ axis, the delay time of the wave arrival at the distant support is $t_0=L$/$v_s=$ 0.86-1.29 s, i.e. $t_0=$ (2.14-3.2)∙$T_{prev}$.

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. 9The design model. The red arrow depicts the direction of a seismic input motion

Fig. 10Accelerogram 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. 11The 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 $t_0=$ 0 and $t_0=$ 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 ($\gamma =\delta$/$\pi =$ 0.1, $\gamma$ is the coefficient of internal inelastic resistance, $\delta$ is the logarithmic decrement of oscillations), and the upper two floors are made of steel ($\gamma =$ 0.03). Soil compliance in the direction of the translational displacements is modeled by the springs with stiffness $k_x=k_y=k_z=$ 10 MN/m, installed in the level of the supports, and the soil dissipation is modeled by local dampers with viscosity ${\beta }_x={\beta }_y={\beta }_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 $f\left(t\right)=1-\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{\Omega }t\right)$, where the frequency of the forced oscillations is $\mathrm{\Omega }={\omega }_2$, and ${\omega }_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:

where ${\mathbf{C}}_e$, ${\mathbf{K}}_e$ are dissipation and stiffness matrices of the finite element $e$, ${\gamma }_e$ is the coefficient of internal inelastic resistance of the material.

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

Fig. 13Horizontal 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.