Numerical Simulations of LEAP Dynamic Centrifuge Model Tests for Response of Liquefiable Sloping Ground

This paper presents numerical simulations (modified Type-B) related to LEAP-UCD-2017 (Liquefaction Experiments and Analysis Projects) dynamic centrifuge model tests for a liquefiable sloping ground conducted by various institutions. The numerical simulations are performed using a pressure-dependent constitutive model implemented with the characteristics of dilatancy and cyclic mobility. The soil parameters are determined based on a series of available stress-controlled cyclic triaxial tests during the Type-A simulation phase for matching the liquefaction strength curves of Ottawa F-65 sand. The computational framework for the dynamic response analysis is discussed. Computed results are presented for the selected centrifuge experiments (modified Type-B simulations). Measured time histories (e.g., displacement, acceleration and excess pore pressure) of these experiments are reasonably captured. Comparisons between the numerical simulations and measured results showed that the pressure-dependent constitutive model as well as the overall employed computational framework have the potential to predict the response of the liquefiable sloping ground, and subsequently realistically evaluate the performance of an equivalent soil system subjected to seismically induced liquefaction.


LEAP (Liquefaction Experiments and Analysis Projects)
is an effort to facilitate validation and verification of numerical procedures for liquefaction-induced lateral spreading analysis of a liquefiable sloping ground (Kutter et al. , 2015Manzari et al. 2014). In order to obtain a set of reliable centrifuge test data with high quality among different centrifuge facilities, a centrifuge experiment was repeated at six centrifuge facilities in LEAP-GWU-2015 (Kutter et al. 2017;Manzari et al. 2018), as one project within LEAP. In addition, the associated numerical simulations conducted by several predictors (Ueda and Iai 2018;Ziotopoulou 2018;Zeghal et al. 2018;Ghofrani and Arduino 2018;Armstrong 2018) were compared with the measured response from the conducted experiments.
As part of the ongoing LEAP, i.e., LEAP-UCD-2017 ), a new set of dynamic centrifuge tests have been performed to simulate the liquefactioninduced lateral spreading phenomenon in a fully saturated sloping ground. A major goal of LEAP-UCD-2017 is to quantify repeatability of the centrifuge test as well as sensitivity of the response with respect to variations of the soil density and base input excitation.
On this basis, this paper presents the results of numerical simulations for the dynamic centrifuge tests and the emphasis is on calibration based on the Type-B simulation phase, and subsequent minor modification based on knowledge of experiment results during the Type-C phase. All the numerical simulations are performed using a pressure-dependent constitutive model (Yang 2000;Elgamal et al. 2003;Parra 1996;Yang and Elgamal 2002) implemented with the characteristics of dilatancy and cyclic mobility. The soil parameters are determined based on a series of stress-controlled cyclic triaxial tests provided during the Type-A simulation phase for matching the liquefaction strength curves of Ottawa F-65 sand. To better capture the overall dynamic response of each selected centrifuge test in the modified Type-B simulations, only one parameter was adjusted by observations from all selected centrifuge test results provided in the Type-C simulation phase (to increase the rate of pore pressure build-up).
The following sections of this paper outline: (1) computational framework, (2) specifics and calibration processes, (3) details of the employed FE modeling techniques, and (4) computed results of the selected centrifuge tests. Finally, a number of conclusions are presented and discussed.

Brief Summary of the Centrifuge Tests
A schematic representation of the centrifuge tests is shown in Fig. 26.1. The soil specimen is a sloping layer of Ottawa F-65 sand with 5 slope (target dry density is 1652 kg/m 3 , D r ¼ 65%). The soil layer has a length of 20 meters (in prototype scale) and a height of 4 m (in prototype scale) at the center. The specimen is built in a container with rigid walls. Three arrays of accelerometers and pore pressure transducers are placed in the central section and at 3.5 m away from the side walls on the right and left of the centrifuge model. In the vertical direction, the sensors are 0.5 m apart. The sensors in the central section are required sensors and those of side arrays (away from central section) are recommended for all centrifuge facilities. Table 26.1 shows a summary of the experiments selected for LEAP-UCD-2017 Type-B simulations including centrifugal acceleration and the achieved density of the soil. All centrifuge models were subjected to a target motion of ramped, 1 Hz sine wave base motion with amplitude 0.15 g. Figure 26.2 shows the target motion and achieved motions for the selected centrifuge experiments.

Finite Element Model
A two-dimensional FE mesh with element size 0.5 m ( Fig. 26.3) is created to represent the centrifuge model with rigid walls. All numerical simulations for the selected nine centrifuge experiments (  (Yang 2000;Yang and Elgamal 2002) under static and seismic loading.
Four-node plane-strain elements with two-phase material following the u-p (Chan 1988) formulation were employed for simulating saturated soil response, where u is the displacement of the soil skeleton and p is the pore water pressure. Implementation of the u-p element is based on the following assumptions: (1) small deformation and rotation; (2) solid and fluid density remain constant in time and space; (3) porosity is locally homogeneous and constant with time; (4) soil grains are incompressible; (5) solid and fluid phases are accelerated equally. Hence, the soil layers represented by effective stress fully coupled u-p elements are capable of accounting for soil deformations and the associated changes in pore water pressure.

Soil Constitutive Model
The employed soil constitutive model (Yang 2000;Elgamal et al. 2003;Parra 1996;Yang and Elgamal 2002) were developed based on the multi-surface-plasticity theory (Mroz 1967;Iwan 1967;Prevost 1978Prevost , 1985. In this employed soil constitutive model, the shear-strain backbone curve was represented by the hyperbolic relationship with the shear strength based on simple shear (reached at an octahedral  shear strain of 10%). The small strain shear modulus under a reference effective confining pressure p 0 r is computed using the equationG ¼ G 0 ( p 0 /p 0 r ) n , where p 0 is effective confining pressure. The dependency of shear modulus on confining pressure is taken as (n ¼ 0.5). The critical state frictional constant M f (failure surface) is related to the friction angle ϕ (Chen and Mizuno 1990) and defined as M f ¼ 6sinϕ/ (3-sinϕ). As such, the soil is simulated by the implemented OpenSees material PressureDependMultiYield02. Brief descriptions of this soil constitutive model are included below.

Yield Function
The yield function is defined as a conical surface in principal stress space (Prevost 1985;Lacy 1986;Yang and Elgamal 2002): where, s ¼ σ 0 À p 0 δ is the deviatoric stress tensor, σ 0 is the effective Cauchy stress tensor, δ is the second-order identity tensor, p 0 is mean effective stress, p 0 0 is a small positive constant (0.3 kPa in this paper) such that the yield surface size remains finite at p 0 ¼ 0 for numerical convenience and to avoid ambiguity in defining the yield surface normal to the yield surface apex. a is a second-order deviatoric tensor defining the yield surface center in deviatoric stress subspace, M defines the yield surface size, and ":" denotes doubly contracted tensor product.

Contractive Phase
Shear-induced contraction occurs inside the phase transformation (PT) surface (η < η PT ), as well as outside (η > η PT ) when _ η < 0, where, η is the stress ratio and η PT is the stress ratio at phase transformation surface. The contraction flow rule is defined as (Yang et al. 2008): where c 1 , c 2 and c 3 are non-negative calibration constants, γ d is octahedral shear strain accumulated during previous dilation phases, p a is atmospheric pressure for normalization purpose, and _ s is the deviatoric stress rate. The _ n and _ s tensors are used to account for general 3D loading scenarios, where, _ n is the outer normal to a surface. The parameter c 3 is used to represent the dependence of pore pressure buildup on initial confinement (i.e., K σ effect).

Dilative Phase
Dilation appears only due to shear loading outside the PT surface (η > η PT with _ η > 0), and is defined as (Yang et al. 2008): where d 1 , d 2 and d 3 are non-negative calibration constants, and γ d is the octahedral shear strain accumulated during all dilation phases in the same direction as long as there is no significant load reversal. It should be mentioned that γ d accumulates even if there are small unload-reload phases, resulting in increasingly stronger dilation tendency and reduced rate of shear strain accumulation.

Neutral Phase
When the stress state approaches the PT surface (η ¼ η PT ) from below, a significant amount of permanent shear strain may accumulate prior to dilation, with minimal changes in shear stress and p 0 (implying p 00 ¼ 0). For simplicity, P 00 ¼ 0 is maintained during this highly yielded phase until a boundary defined in deviatoric strain space is reached, and then dilation begins. This yield domain will enlarge or translate depending on load history (Yang et al. 2008).
It should be noted that PressureDependMultiYield02 has been improved with new flow rules in order to better capture contraction and dilation in sands and the model parameters were calibrated with established guidelines on the liquefaction triggering logic, i.e., cyclic stress ratio versus number of equivalent uniform loading cycles in undrained (direct simple shear) DSS loading to cause single-amplitude shear strain of 3% (Khosravifar et al. 2018).

Boundary and Loading Conditions
The boundary and loading conditions for the dynamic analysis of the sloping ground under an input motion are implemented in a staged fashion as follows: 1. Gravity was applied to activate the initial static state (Fig. 26.4) for the sloping ground with: (a) linear elastic properties (Poisson's ratio of 0.47), (b) nodes on both side boundaries (vertical faces) of the FE model were fixed against longitudinal translation, (c) nodes were fixed along the base against vertical translation, (d) water table was specified ( Fig. 26.2) with related water pressure and nodal forces specified along ground surface nodes. At the end of this step, the static soil state was imposed and displacements under own-weight application were re-set to zero using the OpenSees command InitialStateAnalysis.
2. Soil properties were switched from elastic to plastic. 3. Nodes were fixed along the base against longitudinal translation. 4. Dynamic analysis is conducted by applying an acceleration time history to the base of the FE model.
The FE matrix equation is integrated in time using a single-step predictor multicorrector scheme of the Newmark type (Chan 1988;Parra 1996) with integration parameters γ ¼ 0.6 and β ¼ 0.3025. The equation is solved using the modified Newton-Raphson method, i.e., Krylov subspace acceleration (Carlson and Miller 1998) for each time step. A relatively low level of stiffness proportional damping (coefficient ¼ 0.003) with the main damping emanating from the soil nonlinear shear stress-strain hysteresis response was used to enhance the numerical stability of the liquefiable sloping system. The tolerance criteria used to check the convergence is based on the increment of energy with a tolerance of 10 À6 . The key parameters for calibrating the provided triaxial data are the small strain shear modulus (G o ), the contraction parameters (c 1 , c 2 ), and the dilation parameters (d 1 , d 2 ). For the rest of 16 parameters (e.g., c 3 and d 3 ), default values were used. The internal friction angle and phase transformation angle are determined directly from the laboratory data. The saturated soil mass density is detemined based on the specific gravity of Ottawa F-65 sand and the void ratio. The parmeability of the the Ottawa F-65 sand is 1.1 Â 10 À4 m/s determined from El Ghoraiby et al. (2019). On this basis, the calibrated soil model parameters for the Ottawa F-65 sand at three different void ratios are presented in Table 26.2. Figure 26.5 shows the plot of the liquefaction strength curves for Ottawa F-65 sand with initial void ratios of 0.515, 0.542 and 0.585 (ρ d ¼ 1744.2, 1712.6, 1665.6 kg/m 3 ). The data plotted is composed of the number of cycles until a 2.5% single amplitude of strain is achieved versus the applied cyclic stress ratio. It can be seen that a comparatively good match is reached between simulation results and laboratory test data for the Ottawa sand F-65 at three different void ratios. Figure 26.6 shows an example of comparison results between simulation and experiment for Ottawa F-65 sand with void ratio 0.585, which will be further used in the centrifuge tests of Type-B simulations. It is noted that the yield surface of PressureDependMultiYield02 is conical in principal stress space, resulting in symmetric response for deviatoric stress versus axial strain (Fig. 26.6b). However, it still can reasonably capture the stress path, maximum axial strain, and pore pressure generation (Fig. 26.6a, c, and d).

Computed Results of Type-B Simulations
In this section, the numerical results of selected nine centrifuge tests in the modified Type-B simulations for Ottawa F-65 sand with void ratio e o ¼ 0.585 are presented. The achieved motion for each centrifuge test is shown in Fig. 26.2. As discussed above, the model parameters are calibrated for matching the liquefaction strength curves for the Ottawa F-65 sand with different void ratios from the Type-A simulation phase. Although the measured and simulated results are in good agreement for the case of void ratio e o ¼ 0.585, it was found that the simulation was unable to capture the overall response for the selected nine centrifuge tests in Type-B simulations (once the centrifuge test results were made available for Type-C simulations). In order to better capture the overall dynamic response for all selected centrifuge tests, the contraction parameter c 1 equal to 0.015 in Table. 26.2 was adjusted to c 1 ¼ 0.08 based on the observations from the provided test results in the Type-C simulation phase. Figure 26.7 shows comparison results between numerical simulation and triaxial experiment for Ottawa F-65 sand with void ratio 0.585 (Fig. 26.6) after the contraction parameter c 1 was adjusted. As such, some discrepancies of the soil response can be found, e.g., the pore pressure in Fig. 26.7d generates faster than that in Fig. 26.6d, and the axial strain in Fig. 26.7b is overestimated compared to that in Fig. 26.6b. As discussed above, the contraction parameter c 1 was adjusted to better capture the overall response of the centrifuge tests. Therefore, the following numerical results for each centrifuge test are computed based on the parameters in Table 26.2 with the adjusted c 1 parameter as described above.
The numerical results for UCD-3 are presented in Fig. 26.8. The computed time histories of acceleration (AH1-AH4) and excess pore pressure (P1-P4) reasonably match with those from measurements ( Fig. 26.8a, b). Both the computed results and the measurements showed a consistent trend for negative spikes of pore pressure generation and acceleration due to dilation. The computed lateral displacement in Fig. 26.8c at the midpoint of the ground surface shows a trend similar to the measurement. Figure 26.9 shows the computed horizontal and vertical displacement contours at the end of shaking. It is noted that the arrows in Fig. 26.9a shows the soil movement and the length represents the total displacement. Figure 26.10 shows the computed soil response including the shear stress versus mean effective stress (confinement) and shear stress versus shear strain for the integration points near the locations of pore pressure transducer (P1-P4). It is observed that the pore pressure dips are consistent with the spikes in shear stress due to dilation.
The numerical results for UCD-1 are shown in Fig. 26.11. Comparison of the results in Fig. 26.11a shows that the computed acceleration is in good agreement with that of measurement. In Fig. 26.11b, the time histories of pore pressure are overestimated by numerical simulations. Figure 26.11c shows the comparison of displacement between numerical and measured computed results. It can be seen that the measured displacement is significantly smaller than that of simulation. Comparing the achieved response of UCD-1 with those of other centrifuge tests, the reason for discrepancies between the simulations might stem from the fact that the soil near the ground surface was built denser than what was expected. Figure 26.12 shows the results of NCU-3. The computed time histories of acceleration (AH3-AH4) at deeper depth are in good agreement with those from measurements ( Fig. 26.12a). The acceleration at AH1 and AH2 show very high negative spikes due to dilation, that the computed results did not reproduce. By comparing time histories of pore pressure, a relatively good agreement is shown in Fig. 26.12b in which, time histories of P3 at depth 2 m was not available. Figure 26.12c indicates that the computed lateral displacement at midpoint of the ground surface is accumulating cycle-by-cycle with residual lateral displacement around 0.055 m. Since the measured time history of lateral displacement not being available, the dynamic component computed from the acceleration history at A4 is shown (no residual lateral displacement induced by the inclination of the ground). The results for CU-2 are shown in Fig. 26.13. The simulation results provide a similar response for AH1 and AH2 (Fig. 26.13a). Compared to the time histories of the numerical simulations, the experiment results show shifted negative spikes in acceleration due to dilation. The overall trend of the time histories of the excess pore pressure is generally consistent with the measurements (Fig. 26.13b). In Fig. 26.13c, the computed lateral displacement is around 0.18 m, however, the lateral displacement at the midpoint of the ground in the experiment shows around 0.42 m at the end of shaking. By comparing the displacement with those of the other testing facilities, the induced lateral displacement of CU-2 is relative higher. This might be due to that the soil near the ground surface was built looser than what was expected. Figure 26.14 presents the results of Ehime-2. The computed time histories of acceleration and excess pore pressure are in good agreement with those from measurements ( Fig. 26.14a, b). Although, both the computed results and the measurements showed a consistent trend for the negative spikes of pore pressure generation and acceleration due to dilation, some discrepancies are found in the values of excess pore pressure between the measured and computed results. At the locations of P2-P4, the values of excess pore pressure for the numerical results are slightly lower than those of the measurement, but P1 is slightly higher. The computed lateral displacement in Fig. 26.14c at the midpoint of the ground surface shows a good trend similar to the measurement and with similar residual displacement at the end of shaking.
The results for ZJU-2 are shown in Fig. 26.15. The simulation results provide a similar response for AH1 and AH2 (Fig. 26.15a). The comparison results for ZJU-2 are similar to CU-2 (i.e., the experiment results show shifted negative spikes in acceleration compared to the numerical results). The overall trend of the time histories of excess pore pressure is generally consistent with the measurements (Fig. 26.15b). Discrepancies are also found in the values of excess pore pressure between the measured and computed results. At the locations of P1-P2, the values of excess pore pressure for the numerical results are slightly higher than those of the measurement. In Fig. 26.15c, the computed lateral displacement is around 0.12 m and the lateral displacement at the midpoint of the ground in the experiment shows around 0.3 m at the end of shaking. The numerical results for KAIST-1 are presented in Fig. 26.16. The computed time histories of acceleration and excess pore pressure are in good agreement with those from measurements ( Fig. 26.16a, b). The time histories of AH3 and AH4 in measurement show more negative spikes due to dilation. Both the computed results and the measurements showed a consistent trend for the negative spikes of pore pressure generation due to dilation. In addition, the computed time histories of excess pore pressure at P1 and P2 reasonably captured the pore pressure dispassion at the end of shaking from 18 to 24 s. Figure 26.16c indicates that the computed lateral displacement at the midpoint of the ground surface is accumulating cycle-by-  Fig. 26.16c. Figure 26.17 presents the results of KAIST-2. The computed time histories of acceleration and excess pore pressure are in good agreement with those from measurements ( Fig. 26.17a, b). Compared to the time histories of the numerical simulations, the experiment results show shifted negative spikes in acceleration due to dilation. The overall trend of time histories of excess pore pressure is generally consistent with the measurements (Fig. 26.17b). The computed lateral displacement in Fig. 26.17c at the midpoint of the ground surface shows a good trend similar to the measurement and with similar residual displacement at the end of shaking. Figure 26.18 shows the results of KyU-3. The computed time histories of acceleration (AH1-AH3) are in good agreement with those from measurements ( Fig. 26.18a). The acceleration at AH4 near the ground surface showed very high negative spikes due to dilation and was not captured by the computed result. By comparing the time histories of pore pressure, a relatively good agreement is shown. Figure 26.18c indicates that the computed lateral displacement at the midpoint of the ground surface is accumulating cycle-by-cycle and with residual lateral displacement around 0.055 m. Similar to NCU-3 and KAIST-1, the dynamic component of displacement computed from the acceleration history at AH4 is shown in Fig. 26.18 (no residual lateral displacement).

Conclusions
The numerical simulations (modified Type-B) in LEAP-UCD-2017 centrifuge model tests for a liquefiable sloping ground conducted by various institutions are presented in this paper. The computations are performed using a pressure-dependent constitutive model (PressureDependMultiYield02) implemented with the characteristics of dilatancy and cyclic mobility. The soil parameters are determined based on a series of available stress-controlled cyclic triaxial tests (provided during the Type-A simulation phase) for matching the liquefaction strength curves. The computational framework and staged analysis procedure are presented as well. The primary conclusions can be drawn as follows: 1. Determination of soil model parameters in the Type-A simulation phase showed that the pressure-dependent model, i.e., PressureDependMultiYield02 material, generally resonably captured the liquefaction strength curves for the cyclic stresscontrolled triaxial tests of the Ottawa F-65 sand at three different void ratios. Although the curve of deviatoric stress versus axial strain is symmetric due to the fact that the yield surface defined in principal stress space is conical, the soil model still can resonably capture the stress path, maximum axial strain, and the pore presssure generation for each triaxial cyclic stress-controlled test. With the knowledge of centrifuge test results provided in the Type-C simulation phase, only the contraction parameter c 1 (controlling the pore pressure generation) was adjusted to better capture the dynamic response of the selected centrifuge tests, leading to what is presented herein as the modified Type-B simulations. The adjustment of this parameter is based on observations for all selected centrifuge test results provided in the Type-C simulation phase. This indicates that, to achieve better fits to the experimental data of the liquefiable sloping ground, the soil constitutive model parameters should be calibrated with some additional knowledge from the response of the centrifuge tests, as well as the triaxial tests on the soils. 3. Although the centrifuge tests in various institutions are not providing completely consistent results, measured time histories of acceleration and pore pressure for these experiments are reasonably captured by the numerical simulations, using the same soil constitutive model parameters. In addition, time histories of accumulated displacement are also captured to a reasonable level. Comparisons between these numerical simulations and measured results for each centrifuge test demonstrated that the PressureDependMultiYield02 soil model as well as the overall employed computational methodology have the potential to predict the response of liquefiable sloping ground, and subsequently realistically evaluate the seismic performance of analogous soil systems subjected to seismically induced liquefaction.
Acknowledgments The authors are grateful for the kind invitation by Professors Majid T. Manzari, Mourad Zeghal, and Bruce L. Kutter to participate in LEAP-UCD-2017. Partial funding of this effort was provided by Caltrans through a project in which Ottawa sand is being used for experimentation.