Comparative Studies between Frequency Domain Analysis and Time Domain Analysis on Free Field One-Dimensional Shear Wave Propagation

1. Introduction

In Korea, the Wolsong Low and Intermediate Level Radioactive Waste Disposal Center(WLDC) has been under construction with a total capacity of 800,000 drums. The 1st phase of the construction, which is the underground silo for LILW disposal facilities with 100,000 drum capacity, was completed in 2014 (see Figure 1) [1]. The facility in the 1st phase was constructed underground, 130 meters below sea level consisting of six silos, 23.6 meters in diameter and 50 meters in height (see Figure 2) [2]. Many studies have been conducted and published to verify the safety assessment since the facility was completed [3,4,5].

Large-scale earthquakes occurred in 2016 and 2017, respectively, in Gyeongju and Pohang areas near WLDC, and interest in the stability of the underground silo increased significantly [6,7]. However, unfortunately, there is little research on the seismic analysis of rocks around underground silos for LILW disposal facilities [8,9,10].

In this study, free-field one-dimensional analyses have been performed before the three-dimensional silo dynamic analyses considering the followings:

• Closed-form solutions (SHAKE, SHAKE91, DEEPSOIL, etc) in the frequency domain are available for one-dimensional shear wave propagation in the linearly viscous elastic system subjected to base accelerations.
• Numerical finite element solutions as the time domain analysis can be directly compared to such closed-form solutions in the free-fields including lateral boundary so that we can assess the accuracy of numerical solutions.

In a parallel study to this paper, a new form of finite element equilibrium equation is derived in terms of relative motions, which is essentially the same equation as that expressed in terms of total motions where the base shear force is applied for earthquake load [11]. The accuracy of conventional finite element solutions is evaluated by directly comparing them with closed-form solutions by frequency domain analysis such as SHAKE91 [12].

2. Linear Frequency Domain Analysis

2.1. Main Algorithm of the Frequency Domain Analysis

Frequency domain analysis has been used for the solution of site responses subjected to vertically propagating shear waves as schematically shown in Figure 3. For such analysis, SHAKE [13] has been the most popular computer program because of its simplicity and practicality in using the program. Since SHAKE, more recent versions have been written to improve the user interface and to show graphical outputs such as SHAKE91 [13] and SHAKE2000 [14].

The main characteristics of the wave motions in the horizontally layered system may be described such as in the following statements. In each layer, the horizontal particle motion consists of the upward incident wave and the downward reflected wave. On the interface between the adjacent layers, displacements and stresses are continuous. On the top surface, the amplitude of the incident wave is the same as that of the reflected wave since the shear stresses should be zero on such a free ground surface. Thus, the amplitude on the top surface is equal to twice the magnitude of the incident wave. On the bottom surface, the downward reflected wave is absorbed into the elastic half-space so that the upward incident wave will not be interrupted by the overlying soil deposit. It should be noted that such an upward incident wave is half the magnitude of "outcrop bedrock motion" for the same reason as explained for top ground surface.

The main algorithms of the frequency domain analysis may be described in the following way. For each harmonic motion, set up transfer functions for the incident and reflected waves in each layer, refer to SHAKE for the detailed derivation. These transfer functions represent the ratio of amplitudes in a layer to those at the top surface. The input object accelerations in the time domain are converted to Fourier series form in the frequency domain using the Fast Fourier Transform (FFT) method. Amplitudes at any location in the layer can be found by using Fourier series and transfer functions in the frequency domain and then responses in the time domain can be determined by inversing FFT.

2.2. Site Profile and Response Spectra for Input Earthquake

Figure 4 and Table 1 show the site profile and material properties used for all free-field analyses for both frequency domain and time domain analyses.

Figure 5 shows input acceleration time history which has been used for the three-dimensional silo dynamic analyses as well as one-dimensional free-field analyses. This is the earthquake of local magnitude of ML = 5.8 recorded at the MKL station in Gyeongju Myeonggye-ri, South Korea on September 12, 2016 [7]. This input acceleration is applied as outcrop motion at the base rock located 373 meters below the ground surface.

The characteristics of this earthquake can be viewed by plotting the response spectra as shown in Figure 6. It shows that absolute spectral accelerations are so high in periods between 0.03 and 0.06 seconds, indicating peak spectra values occurring at such high-frequency regions.

Two computer programs (SHAKE91 [13] and DEEPSOIL [15]) are selected for the closed-form solutions in the frequency domain. Figure 7 shows the portion of the input user interface captured from the program DEEPSOIL.

2.3. Analysis Results

Figure 8 shows the acceleration time histories on the ground surface which compare SHAKE91 to DEEPSOIL. Figure 9 shows the same comparison between t = 2 and t = 4 seconds, showing almost identical responses in the time where strong motions occur.

Figure 10 shows the acceleration time histories at a depth of 138 meters which compare SHAKE91 to DEEPSOIL. Figure 11 shows the same comparison between t = 2 and t = 4 seconds, showing almost identical responses in the time when strong motions occur. Note that this depth corresponds to the location of the silo mid height as shown in Figure 4.

3. Linear Time Domain Analysis

3.1. Finite Element Formulations of Dynamic Equation

In the previous section, the frequency domain analysis solves the site responses based on wave equations with appropriate boundary conditions using the Fourier transformation method which converts input motions from the time domain to the frequency domain and then gets responses in the time domain by inverse transformation.

Time domain analysis involves setting up the global dynamic equilibrium in a spatially discretized system and solving the equilibrium equation directly by the time marching schemes such as Newmark method. Finite element analysis is the most popular and powerful method that has been used to obtain dynamic solutions in the time domain.

Finite element formulations of dynamic equations for linearly viscous elastic material subjected to base earthquake motions may be expressed in the following matrix form at time step n.

$$M {\ddot{u}}_n+D {\dot{u}}_n+K u_n=R_n$$

(1)

$$R_n=-M · I · {\ddot{u}}_{gn}\mathrm{when}u\mathrm{represents}\mathrm{total}\mathrm{displacement}$$

$$R_{\mathrm{n}}={\mathsf{\rho }}_{\mathrm{r}}·{\mathrm{c}}_{\mathrm{sr}}·{\mathrm{A}}_{\mathrm{s}}·J · {\dot{\mathrm{u}}}_{\mathrm{gn}}\mathrm{when}u\mathrm{represents}\mathrm{relative}\mathrm{displacement}$$

$$I{}^{\mathrm{T}}=〈1 1 ··11〉$$

$$J{}^{\mathrm{T}}=〈0 0 ··01〉$$

Where

• $M, D, K :\mathrm{Mass},\mathrm{damping}\mathrm{and}\mathrm{stiffness}\mathrm{matrices}$
• ${\dot{\mathrm{u}}}_{\mathrm{gn}},{\ddot{\mathrm{u}}}_{\mathrm{gn}} : \mathrm{Earthquake}\mathrm{outcrop}\mathrm{velocity}\mathrm{and}\mathrm{acceleration}\mathrm{at}\mathrm{bedrock}$
  • ${\mathsf{\rho }}_{\mathrm{r}}, {\mathrm{c}}_{\mathrm{sr}}:\mathrm{Mass}\mathrm{density}\mathrm{and}\mathrm{shear}\mathrm{wave}\mathrm{velocity}\mathrm{at}\mathrm{bedrock}$
  • ${\mathrm{A}}_{\mathrm{s}} : \mathrm{Tributary}\mathrm{area}\mathrm{where}\mathrm{base}\mathrm{shear}\mathrm{force}\mathrm{is}\mathrm{acting}$

For the direct time integration of Eq. (1), Newmark constant average acceleration method, which is unconditionally stable, may be used.

$${\dot{u}}_{\mathrm{n}}={\dot{u}}_{\mathrm{n}-1}+\left({\ddot{u}}_{\mathrm{n}-1}+{\ddot{u}}_{\mathrm{n}}\right)·\left(\mathsf{\Delta }\mathrm{t}/2\right)$$

(2)

$$u_{\mathrm{n}}=u_{\mathrm{n}-1}+{\dot{u}}_{\mathrm{n}-1}·\mathsf{\Delta }\mathrm{t}+\left({\ddot{u}}_{\mathrm{n}-1}+{\ddot{u}}_{\mathrm{n}}\right)·\left({\mathsf{\Delta }\mathrm{t}}^2/4\right)$$

(3)

From Eqs. (2) and (3), we can obtain the following equations.

$${\ddot{u}}_{\mathrm{n}}=\left(4/{\mathsf{\Delta }\mathrm{t}}^2\right)·u_{\mathrm{n}}-{\mathrm{A}}_{\mathrm{n}-1}$$

(4)

$${\dot{u}}_{\mathrm{n}}=\left(2/\mathsf{\Delta }\mathrm{t}\right)·u_{\mathrm{n}}-{\mathrm{B}}_{\mathrm{n}-1}$$

(5)

Substituting Eqs. (4) and (5) into Eq. (1), we obtain the following equation.

$$\overline{K} ·u_{\mathrm{n}}={\overline{R}}_{\mathrm{n}}$$

(6)

$$\overline{K}=\left(4/{\mathsf{\Delta }\mathrm{t}}^2\right)·M+\left(2/\mathsf{\Delta }\mathrm{t}\right)·D+K$$

(7)

$${\overline{R}}_{\mathrm{n}}= R_{\mathrm{n}}+{\mathrm{A}}_{\mathrm{n}-1}·M+{\mathrm{B}}_{\mathrm{n}-1}·D$$

(8)

$${\mathrm{A}}_{\mathrm{n}-1}=4·\left(u_{\mathrm{n}-1}/{\mathsf{\Delta }\mathrm{t}}^2+{\dot{u}}_{\mathrm{n}-1}/\mathsf{\Delta }\mathrm{t}+{\ddot{u}}_{\mathrm{n}-1}/4\right)$$

(9)

$${\mathrm{B}}_{\mathrm{n}-1}=2·\left(u_{\mathrm{n}-1}/\mathsf{\Delta }\mathrm{t}\right)+{\dot{u}}_{\mathrm{n}-1}$$

(10)

It is quite common practice to express the viscous damping by Rayleigh and Lindsay [16] which consists of mass and stiffness proportional terms at the element level.

$$D=\mathrm{a}·M+\mathrm{b}·K$$

(11)

The energy loss associated with such viscous damping in Eq. (11) is proportional to the velocity, which is also dependent on the frequency of the motion. The energy dissipation in soils, however, is independent of frequency even at very small strain levels based on experimental test data [17].

To mitigate such a frequency dependency in Rayleigh damping, the values of "a" and "b" in Eq. (11) are expressed in terms of two target frequencies $\left({\mathsf{\omega }}_1{\mathrm{and}\mathsf{\omega }}_{\mathrm{i}}\right)$ [18].

$$\mathrm{a}=2·\mathsf{\beta }·{\mathsf{\omega }}_1·{\mathsf{\omega }}_{\mathrm{i}}/\left({\mathsf{\omega }}_1+{\mathsf{\omega }}_{\mathrm{i}}\right)$$

(12)

$$\mathrm{b}=2·\mathsf{\beta }/\left({\mathsf{\omega }}_1+{\mathsf{\omega }}_{\mathrm{i}}\right)$$

(13)

where

${\mathsf{\omega }}_1:$	Fundamental natural circular frequency of the system
${\mathsf{\omega }}_{\mathrm{i}}:$	Predominant circular frequency of the input earthquake motion
$\mathsf{\beta }:$	Critical damping ratio in an element

3.2. Computer Programs for Free Field Analysis

Three finite element computer programs are selected to perform the free-field analysis in the time domain; SRAP-1D, QUAD-4M, and SMAP-3D. SRAP-1D is a one-dimensional program that is specific to site response analysis subjected to earthquake motions. The original lumped mass solution listed in Ref. [19] is modified to include elastic half-space. QUAD-4M is a two-dimensional program to evaluate the seismic response of soil structures [20]. SMAP-3D is the general-purpose three-dimensional program developed by Comte Research [21].

All three programs use the following schemes for free-field analyses:

• Lumped mass matrix as implicitly used by the frequency domain analysis
• Transmitting boundary on bottom of finite element mesh to represent the elastic half-space
• Two methods for applying external earthquake loadings [11]
  • Method 1: Base accelerations to relative displacement fields as the conventional procedure
  • Method 2: Base shear forces associated with base velocities
• Newmark average acceleration method for the time integration

Figure 12 shows the cross-section of the 3D finite element mesh modeled by SMAP-3D.

3.3. Analysis Results

Figure 13 shows the acceleration time histories on the ground surface which compare the results of all three computer programs. Figure 14 shows the same comparison between t = 6 and t = 8 seconds, showing almost identical responses in the time where strong motions occur. It also shows that Method 1 conventional procedure produces almost the same results as Method 2. It should be noted that Method 1 would produce exactly the same results as Method 2.

Figure 15 shows the acceleration time histories at a depth of 138 meters which compare QUAD-4M to SMAP-3D. Figure 16 shows the same comparison between t = 6 and t = 8 seconds, showing almost identical responses in the time when strong motions occur. Note that this depth corresponds to the location of the silo mid height as shown in Figure 4.

4. Linear Frequency vs Linear Time Domain Analysis

4.1. General

In this section, SMAP-3D time domain solutions are compared to the exact closed-form frequency domain SHAKE91 solutions for the analysis of free-field responses subjected to vertically propagating shear waves caused by earthquake motions.

As described in the previous Section 3, finite element time domain solutions employ the frequency-dependent Rayleigh formulation to model material damping. Consequently, responses would be more damped at some frequencies and less damped at the other frequencies. On the other hand, the damping scheme in SHAKE91 is independent of frequency.

Both SHAKE91 and SMAP-3D solutions are presented separately in the previous Section 2 and Section 3, respectively. In this section, one additional analysis is conducted for SMAP-3D using numerical damping in Newmark β method with γ = 1. Obviously, such numerical damping will damp down the responses associated with high-frequency modes. However, it is often applied to obtain stable solutions when the material state changes abruptly due to brittle failure or joint separation during strong shaking dynamic motions.

4.2. Analysis Results

Figure 17 shows the acceleration time histories on the ground surface which compare the SMAP-3D with γ = 0.5 to SHAKE91. Figure 18 shows the same comparison between t = 6 and t = 8 seconds. Compared to SHAKE91, SMAP-3D results show somewhat higher responses at the time when strong motions occur.

Figure 19 shows the acceleration time histories at a depth of 138 meters which compares the SMAP-3D with γ = 0.5 to SHAKE91. Figure 20 shows the same comparison between t = 6 and t = 8 seconds, where strong motions occur. SMAP-3D results at this depth show closer to SHAKE91 results than we see those comparisons on the ground surface. Note that this depth corresponds to the location of the silo mid height as shown in Figure 4.

Figure 21 shows the acceleration time histories on the ground surface which compare the SMAP-3D with γ = 1 to SHAKE91. Figure 22 shows the same comparison between t = 6 and t = 8 seconds. Compared to SHAKE91, SMAP-3D results show significantly damped responses in the time when strong motions occur.

Figure 23 shows the acceleration time histories at a depth of 138 meters which compare the SMAP-3D with γ = 1 to SHAKE91. Figure 24 shows the same comparison between t = 6 and t = 8 seconds. Compared to SHAKE91, SMAP-3D results show moderately damped responses in the time when strong motions occur. Note that this depth corresponds to the location of the silo mid height as shown in Figure 4.

So far, all comparisons are based on the acceleration time histories. Now we want to make some comparisons based on relative displacement and shear stress time histories.

Figure 25 shows the relative displacement time histories on the ground surface which compare the SMAP-3D with γ = 0.5 and γ = 1 to SHAKE91. Figure 26 shows the same comparison between t = 7 and t = 9 seconds. Compared to SHAKE91, SMAP-3D results are closer to SHAKE91 results than we see those comparisons for acceleration histories. It is also noticed that the effect of numerical damping is less significant in the relative displacements.

Figure 27 and Figure 28 show the relative displacement time histories at the depth 138 m and on the bottom surface, respectively, which compare the SMAP-3D with γ = 0.5 and γ = 1 to SHAKE91. We can see the same trends as in the ground surface. Here also, the numerical damping has a slight influence on the response of relative displacements.

Figure 29 shows the shear stress time histories at a depth of 138 m which compares the SMAP-3D with γ = 1 to SHAKE91. Figure 30 shows the same comparison between t = 6 and t = 8 seconds. SMAP-3D calculation with numerical damping of γ = 1 predicts reasonably well SHAKE91 closed-form results. Note that this depth corresponds to the location of the silo mid height as shown in Figure 4.

5. Conclusions

Comparative studies between frequency domain analysis and time domain analysis on free field one-dimensional analysis are presented in this paper. Free-field one-dimensional analyses have been carried out considering the following:

• Closed-form solutions (SHAKE, SHAKE91, DEEPSOIL, etc) in the frequency domain are available for one-dimensional shear wave propagation in the linearly viscous elastic system subjected to base accelerations.
• Numerical finite element solutions as the time domain analysis can be directly compared to such closed-form solutions in the free-fields including lateral boundary so that we can assess the accuracy of numerical solutions.

The following conclusions arise from numerical studies presented in this paper.

(1)

Two computer programs (SHAKE91 and DEEPSOIL) are selected for the closed-form solutions in frequency domain. The acceleration time histories on the ground surface which compare SHAKE91 to DEEPSOIL show almost identical responses in the time where strong motions occur. The acceleration time histories at the location of the silo mid-height show also almost identical responses in the time where strong motions occur.

(2)

Three finite element computer programs (SRAP-1D, QUAD-4M, and SMAP-3D) are selected to perform the free-field analysis in the time domain. The acceleration time histories on the ground surface which compare the results of all three computer programs show almost identical responses in the time where strong motions occur. It also shows that Method 1 conventional procedure produces almost the same results as Method 2. The acceleration time histories at the location of the silo mid-height which compare QUAD-4M to SMAP-3D show almost identical responses in the time where strong motions occur.

(3)

SMAP-3D time domain solutions are compared to the exact closed-form frequency domain SHAKE91 solutions for the analysis of free-field responses subjected to vertically propagating shear waves caused by earthquake motions. As a result of comparing SMAP-3D with γ = 0.5 and SHAKE91 to the acceleration time histories on the ground surface, SMAP-3D results show somewhat higher responses in the time where strong motions occur. SMAP-3D results at the location of the silo mid-height show more close to SHAKE91 results than we see those comparisons on the ground surface.

(4)

Compared to SHAKE91, SMAP-3D with γ = 1.0 results show significantly damped responses in the time where strong motions occur to the acceleration time histories on the ground surface. Compared to SHAKE91, SMAP-3D results show also moderately damped responses in the time where strong motions occur to the acceleration time histories at the location of the silo mid-height.

(5)

Compared to SHAKE91, SMAP-3D with γ = 0.5 and γ = 1 results show closer to SHAKE91 results than we see those comparisons for acceleration histories to the relative displacement time histories on the ground surface. It is also noticed that the effect of numerical damping is less significant in the relative displacements. Compared to SHAKE91, SMAP-3D with γ = 0.5 and γ = 1 results in the relative displacement time histories at the location of the silo mid-height showing the same trends as in the ground surface. Here also, the numerical damping has a slight influence on the response of relative displacements.

(6)

SMAP-3D calculation with numerical damping of γ = 1 to the shear stress time histories at the location of the silo mid-height predicts reasonably well SHAKE91 closed-form results.