The Importance of Accurate Time-integration in the Numerical Modelling of P-wave Propagation

The numerical dissipation characteristics of the Newmark and generalised-a time-integration schemes are investigated for P-wave propagation in a fully saturated level-ground sand deposit, where higher frequencies than those for S-waves are of concern. The study focuses on resonance, which has been shown to be of utmost importance for triggering liquefaction due to P-waves alone. The generalised-a scheme performs well, provided that the time-step has been carefully selected. Conversely, the dissipative Newmark method can excessively damp the response, changing radically the computed results. This implies that a computationally prohibiting small time-step would be required for Newmark to provide an accurate solution. In finite element (FE) analysis of wave propagation problems the response needs to be computed in both space and time. For the latter, the use of direct numerical methods has been widely employed to solve the dynamic equilibrium equation. According to Hilber and Hughes [8], an essential characteristic of a time-marching scheme is the possession of numerical damping at higher frequencies. This is necessary to filter spurious spikes resulting from poor spatial discretisation, which would otherwise render the interpretation of the high-frequency response of the system problematic [14]. The most widely used time-integration schemes in the field of geotechnics are the Newmark's scheme [17], the Wilson-h method [25] and the HHT method of the family of a time-marching schemes [9]. Due to its superior ability to control dissipation and second-order accuracy, another version of the 'a' schemes, the generalised-a method of Chung and Hulbert [2], subsequently denoted as CH, which has been widely used in one-phase structural dynamic FE analysis, was modified for the two-phase behaviour of soils and implemented into the Imperial College Finite Element Program (ICFEP, [19]) by Kontoe [13] and Kontoe et al. [14]. Over the last few years, the CH method has been used more extensively in the field of geotechnics to integrate both one-phase and two-phase problems [16,11,20,18,22,15], as well as problems involving partially saturated conditions [3]. The unconditional stability of the CH method for the two phase formulation was further investigated by Han et al. [6]. He also showed that the conditions for unconditional stability are applicable to most of the commonly used schemes, including the Newmark algorithms. Although the performance of the various algorithms is understood for linear free-vibration problems, their performance under more complex forced vibration regimes is not as clear, due to the complexity of the analytical solutions. An extensive numerical …


Introduction
In finite element (FE) analysis of wave propagation problems the response needs to be computed in both space and time. For the latter, the use of direct numerical methods has been widely employed to solve the dynamic equilibrium equation. According to Hilber and Hughes [8], an essential characteristic of a timemarching scheme is the possession of numerical damping at higher frequencies. This is necessary to filter spurious spikes resulting from poor spatial discretisation, which would otherwise render the interpretation of the high-frequency response of the system problematic [14].
The most widely used time-integration schemes in the field of geotechnics are the Newmark's scheme [17], the Wilson-h method [25] and the HHT method of the family of a time-marching schemes [9]. Due to its superior ability to control dissipation and second-order accuracy, another version of the 'a' schemes, the generalised-a method of Chung and Hulbert [2], subsequently denoted as CH, which has been widely used in one-phase structural dynamic FE analysis, was modified for the two-phase behaviour of soils and implemented into the Imperial College Finite Element Program (ICFEP, [19]) by Kontoe [13] and Kontoe et al. [14]. Over the last few years, the CH method has been used more extensively in the field of geotechnics to integrate both onephase and two-phase problems [16,11,20,18,22,15], as well as problems involving partially saturated conditions [3]. The unconditional stability of the CH method for the two phase formulation was further investigated by Han et al. [6]. He also showed that the conditions for unconditional stability are applicable to most of the commonly used schemes, including the Newmark algorithms.
Although the performance of the various algorithms is understood for linear free-vibration problems, their performance under more complex forced vibration regimes is not as clear, due to the complexity of the analytical solutions. An extensive numerical comparison of the performance of the various schemes under such complicated conditions was carried out by Kontoe [13] and Kontoe et al. [14], showing the superiority of the CH time-integration scheme in non-linear elasto-plastic transient problems over other common schemes. The main results indicated that even under such conditions, the CH scheme retains the ability to control algorithmic dissipation of the higher frequency modes without altering the lower frequencies. Additionally, it was found not to be affected by the predominant frequency of the input excitation and the natural frequency of the numerical model, although extra care in the selection of the time-step is required under resonant conditions. Conversely, the widely used dissipative scheme of the Newmark family (b dyn = 0.3025, c dyn = 0. 6 It is interesting to note that all investigations of the performance of the various time-integrations schemes in problems involving earthquake loading have been carried out under the shearing mode. No guidance exists on the use of time-integration schemes when modelling the vertical component of the seismic ground motion, where typically higher frequencies are of concern. Given the above, the current study investigates the performance of the CH and the Newmark methods in problems of vertical propagation of P-waves in site response analyses under resonance conditions. Recently, Tsaparli et al. [24] showed that the vertical component of the ground motion can induce cyclic mobility when resonance takes place. Therefore, the examined problem is a challenging test for the time-integration schemes, as the response involves plasticity, hydro-mechanical coupling and high frequencies. The results provide further evidence of the ability of the CH scheme to retain its advanced characteristics under such conditions, whereas it is shown that NMK2 significantly underperforms and can lead to erroneous conclusions.

Description of numerical analyses
In Tsaparli et al. [24] it was shown that in the case of the vertical seismic component and liquefaction triggering in fully saturated level-ground sand deposits, the occurrence of resonance is of utmost importance, as once triggered, significant amplification of the input motion takes place towards the ground surface. This in turn can result in the development of substantial deviatoric stresses, plasticity and possibly liquefaction. For the purposes of the current study the same geometry, sand properties and input acceleration as in Tsaparli et al. [24] are utilised, with the emphasis placed on the behaviour of two time-integration schemes, the CH scheme with numerical parameters equal to b dyn = 0.3025, c dyn = 0.6, a m = 0.35 and a f = 0.45 and the NMK2 scheme. The selection of the above parameters was done so that both schemes result in a spectral radius in the high frequency range, q 1 , equal to 0.818.

Modelled sand deposit
A hypothetical fully-saturated 40 m deep Fraser River Sand (FRS) deposit with the water table at ground level is used in the simulations. The sand deposit has a relative density of 40%, while the small strain shear modulus follows a non-linear distribution with depth according to the Hardin [7] expression. The assumed material properties for the FRS are listed in Table 1. As the Pwave velocity, V P , is governed by the bulk modulus of the pore water, K w ¼ 2:2 Â 10 6 kPa (Eq. (1)), its average value corresponds to 1611 m/s. This, based on the small strain elastic properties, results in a fundamental frequency for the compressional mode of 10.3 Hz (Eq. (2)). spectra radius of the time-integration scheme at infinity where M and q b are the constrained modulus and bulk density of the soil, respectively, H is the depth of the sand column and N is the vibration mode.

Input ground motion
The sand deposit is subjected to the vertical component (V) of the outcrop M w 6.2 22nd February 2011 Christchurch motion [5], as recorded at the Lyttelton Port Company (LPCC) strong ground motion station (SMS), with a modelled duration of 20 s (Fig. 1 input motion denoted as CV). The motion was chosen as its frequency content is very wide, with significant Fourier amplitudes present close to the fundamental frequency of the 40 m deep deposit for P-waves, i.e. at 10.3 Hz.

Numerical method
Effective stress-based finite element (FE) analyses were carried out with the Imperial College Finite Element Program (ICFEP, [19]). The coupling between the solid skeleton and the pore fluid is modelled using the u-p hydro-mechanical formulation [27]. In Tsaparli et al. [24] it was noted that for the sand permeability value and characteristics of the ground motion and deposit used in the study the problem under investigation is within the range over which the u-p formulation is valid [27]. Plane strain conditions are modelled, with the mesh consisting of a single column of 160 Â 1 8-noded quadrilateral elements with pore water pressure degrees of freedom at the 4 corner nodes. To satisfy Bathe's [1] recommendations for frequencies up to about 30 Hz and 8-noded solid elements, an element size of 0.25 Â 0.25 m 2 was adopted to ensure that waves of short wavelengths are not filtered out. Tied degrees of freedom are employed at the lateral boundaries to ensure 1D soil response [26], while the horizontal displacements are prescribed to be zero along the bottom boundary of the mesh. The hydraulic regime in the soil column is defined by restricting the flow at the base of the mesh, imposing zero pore water pressure change at the top boundary and tying the pore water pressure degrees of freedom at the lateral nodes to ensure 1D flow and drainage through the ground surface. Finally, the acceleration time-history is applied incrementally in the vertical direction at the bottom boundary. In all analyses the non-linear solver is based on a modified Newton-Raphson scheme [19].
As both schemes are unconditionally stable, the time-step was selected solely based on accuracy considerations, as the smaller value of one tenth of the following two criteria: (i) the natural period of the deposit, i.e. Dt = 0.0097 s, (ii) the period corresponding to the maximum significant frequency of the input excitation expected to excite higher modes of the system, as obtained from the Fourier Spectrum shown in Fig. 1 (e.g. f = 30 Hz), i.e. Dt = 0.0033 s.
The accuracy of a time-integration method is normally quantified in terms of period elongation and amplitude decay (e.g. [1]). The time-step selection herein aimed to result in practically zero period elongation error for the two dissipative schemes, which can successfully filter the spurious high frequency modes as it will be shown in the subsequent section, so that the impact of the amplitude decay error on the computed response is isolated [13,14]. As such, a time-step as small as 0.003 s was chosen. This value was used in all analyses, so as to isolate the numerical performance of the CH and the NMK2 schemes.
Sand behaviour and liquefaction was modelled using a twosurface bounding surface plasticity model, implemented in ICFEP in three-dimensional stress space [21,23]. Details can be found in Taborda [21] and Taborda et al. [23]. A table with the model parameters for Fraser River Sand, as established by Klokidi [12], can be found in Appendix A.

Ground surface acceleration time-histories and Fourier spectra
As the aim is to compare the performance of the algorithms in the high frequency range, the acceleration response is examined. To this end, Fig. 2a presents the ground surface acceleration time-histories for the two different time-integration schemes (i.e. CH and NMK2) together with the input acceleration, while Fig. 2b depicts the ground surface Fourier spectra for the two analyses. The surface Fourier spectrum of an additional analysis using the average acceleration method (b dyn = 0.25, c dyn = 0.5) of the Newmark family (denoted as NMK1) is also included in Fig. 2c. NMK1 is characterised by a spectral radius, q 1 , of unity which corresponds to zero numerical damping and therefore this method usually suffers from spurious oscillations in the high frequency range. The response predicted by NMK1 can be used as a benchmark though, if the high-frequency noise is isolated, as it is considered the most accurate unconditionally stable scheme [4], featuring zero amplitude decay error.
It becomes immediately evident that, for both schemes (i.e. NMK2 and CH), resonance occurs due to the presence of substantial The results also show that NMK2 filters the spurious spikes at high frequencies that are present in the NMK1 results (i.e. the spike at 55 Hz in Fig. 2c), however at the expense of damping the response excessively at the fundamental frequency, thus resulting in much lower amplitudes of acceleration towards the ground surface, as compared to the CH scheme (note that the NMK1 time domain response is not included in Fig. 2a for clarity, because it is completely dominated by the spurious high frequency response). Conversely, the amplification of the CH analysis at the fundamental frequency, f 0 , agrees with that of NMK1, showing its superiority over the NMK2. In particular, the relative amplitude error at f 0 , as defined by Eq. (3), is À70% for the NKM2, while practically zero for the CH scheme. It should be noted that a negative value implies amplitude decay.
where FA n is the Fourier amplitude for any time-integration scheme and FA NMK1 is the Fourier amplitude for the benchmark case. The above findings can be better understood if the spectral radius, q, and amplitude decay, AD, of each time-integration method are plotted against the ratio of time-step over the period of a single degree of freedom system (SDOF), Dt=T (see Fig. 3). These plots are based on the analysis of Hughes [10], among others, for the free vibration of a linear elastic SDOF and can be qualitatively used to interpret the numerical performance of the two integration schemes. The adopted time-step of 0.003 s corresponds to a ratio of 0.03 and 0.09, considering the fundamental period of the deposit and the period corresponding to a frequency of 30 Hz, the maximum significant frequency obtained from the frequency content of the input motion. Based on the spectral radius plot in Fig. 3a some numerical damping appears to have developed in the case of NMK2 for the critical ratio of 0.09. This, however, appears to be much more significant in the AD error plot shown in Fig. 3b, tying in with the excessively damped response observed previously. Conversely, the CH method, due to its second-order accuracy and controllable numerical dissipation of the high frequencies without impacting the lower ones [14], results in a much smoother transition of the spectral radius from a value of 1 to a value of 0.818 at infinity and, hence, to almost no amplitude decay at a Dt=T ratio of 0.09. The above imply that an accurate solution with the NMK2 method would require an even smaller time-step than 0.003 s, resulting in a prohibitive increase of the computational cost: analyses with different time steps were conducted and showed that only when a Dt value an order of magnitude smaller than the original one was used (i.e. 0.0003 s), the analysis with the NMK2 scheme captured the true response. This agrees with the time-step value obtained from Fig. 3b in order to achieve a very small Dt/T value (%0.01) and, hence, practically zero amplitude decay for NMK2 and for a period corresponding to a frequency of 30 Hz.

Stress paths
The impact of the above observations on the analyses' results are further investigated here by looking at the effective stress paths in mean effective -deviatoric stress space (p 0 À J) at depths of 10 and 20 m for the two analyses under consideration (Fig. 4). The effect of the difference in the numerical damping of the two schemes on the observed response is remarkable: in the case of NMK2 the changes in the direct effective stresses are small due  to the previously seen excessive damping at the fundamental frequency of the deposit (Fig. 2b). This results in a practically linear elastic response. On the other hand, in the case of the analysis with the CH scheme, the accuracy of the numerical solution allows significant changes in the direct effective stresses to take place as a result of resonance, leading to the development of substantial deviatoric stresses, plasticity and, ultimately, liquefaction.

Conclusions
The present numerical study focuses on a comparison of the numerical dissipation characteristics of two time-integration schemes in vertical site response analyses. These are a widely used version of the Newmark family of algorithms which possesses numerical damping [17] and the CH generalised a-method [2,14], which has been used to a lesser extent in geotechnics, where problems involving two phase media are common. In particular, the numerical parameters were chosen such that both methods result in a spectral radius at infinity, q 1 , of 0.818. The time-step was also kept the same, so that the effects of the frequency content of the excitation and of the fundamental frequency of the deposit under resonant conditions on the algorithmic dissipation performance can be isolated and examined. This is of interest as no information is currently available on the performance of the above schemes in coupled hydro-mechanical modelling of problems involving P-wave propagation. The results of the present study show that the CH scheme retains its properties of controllable numerical dissipation even at resonant conditions in the high frequency range. Conversely, the Newmark scheme results in excessive damping of the response at the fundamental frequency, meaning that a prohibitive small time-step would be required to achieve an accurate solution. The numerical error in terms of amplitude decay in this latter case can be so large that it can completely prevent the occurrence of plasticity and ensuing liquefaction in a sand deposit as a result of P-wave propagation alone.