LEAP-UCD-2017 Numerical Simulation at Meisosha Corp

As a part of the LEAP-UCD-2017 project, type B simulations were performed to predict the results of nine selected centrifuge model tests using a two-dimensional effective stress analysis program incorporating the cocktail glass model. In Phase I of the numerical simulation, three parameter sets of the cocktail glass model were determined for three different dry densities. In Phase II, numerical models of the ﬁ nite element method were made based on the data of the centrifuge model. At that time, for each centrifuge model test, we decided which parameter set of the cocktail glass model to apply by referring to the dry density of the analyzed ground. There was a tendency that the response values of the numerical simulation were slightly larger than the response values expected from the difference between the achieved dry density of the centrifuge model and the dry density on which the applied parameter set depends.


Introduction
It is considered advantageous that the constitutive model of sand for numerical simulation is based on the interaction between the sand particles. The strain space multiple mechanism model developed by Iai (Iai et al. 2011), which originated from the multi-inelastic spring model proposed by Towhata and Ishihara (1985), is such a model based on the soil particle interaction (Iai and Ueda 2016). The strain space multiple mechanism model has excellent extensibility. The model can be expanded to a three-dimensional constitutive model (Iai and Ozutsumi 2005). The model also can be expanded to a constitutive model for clay (Iai et al. 2015).
We participated as one of the numerical simulation teams of the LEAP-UCD-2017 project and got the opportunity to perform type B simulation targeted at the centrifuge model tests of liquefiable ground with a gently inclined surface. In the numerical simulations, we applied "the constitutive cocktail glass model" (Iai et al. 2011) to reproduce the behavior of sand. The cocktail glass model is one of the achievements of the strain space multiple mechanism model. The cocktail glass model was incorporated in a two-dimensional dynamic effective stress analysis program, "FLIP ROSE," based on the finite element method. The program is able to analyze earthquake responses accompanied by liquefaction by considering the permeability of the ground. We used this program to perform type B LEAP simulation to confirm the performance of this constitutive model and program.

Model Description
The "cocktail glass model" was proposed for sandy soil based on the framework of a strain space multiple mechanism model for granular materials (Iai et al. 2011;Iai et al. 2013).
The macroscopic effective stress in granular materials is given by the average over contact forces P within a representative volume element with volume V as follows (Iai et al. 2013). The framework of the strain space multiple mechanism model is used to upscale the micromechanical structure into the macroscopic structure (Iai et al. 2013).
t n h i¼ cos ω sin ω sin ω À cos ω ! ð25:3Þ Here, p represents isotropic stress and q virtual simple shear stress (Iai et al. 2013). In order to define the relationship between macroscopic strain tensor ε and macroscopic effective stress tensor σ 0 through the structure defined by Eq. (25.2), let volumetric strain ε and virtual simple shear strain γ(ω) be defined as follows (Iai et al. 2013).
where the double dot symbol denotes a double contraction. To account for the effect of the volumetric strain due to dilatancy (ε d ), effective volumetric strain ε 0 is introduced as follows (Iai 1993).
Effective volumetric strain ε 0 and virtual simple shear strain γ are used respectively to define isotropic stress p and virtual simple shear stress q as follows (Iai et al. 2013). Fig. 25.1 Contact between two particles: normal n, tangential direction t, and contact force P (Iai et al. 2013) q Dilatancy ε d is decomposed into contractive dilatancy ε c d and dilative dilatancy ε d d as follows (Iai et al. 2011).
The rate of ε c d is controlled by the irreversible (plastic) portion of virtual simple shear strain rate ( _ γ p ) and the rate of ε d d is controlled by virtual simple shear strain rate ( _ γ ).

Overview
We used a two-dimensional dynamic effective stress analysis program, "FLIP ROSE," in the type B predictions. "FLIP ROSE" (Finite element analysis program of Liquefaction Process/Response of Soil structure systems during Earthquakes) is based on the finite element method and infinitesimal strain formulation. The program provides a multistage analysis function that enables the analysis to reproduce the initial stress state accurately by considering the construction process.

Elements
The program has 20 kinds of element. We used the pore water elements (partially drained) for modeling the behavior of the pore water and the cocktail glass model elements for modeling the behavior of the sand soil skeleton. The pore water element is based on u-p formulation by Zienkiewicz and Bettss (1982). The cocktail glass model element is based on the model proposed by Iai et al. (2011). For simulation of soil behavior below the ground water table, the pore water element (partially drained) is superposed with the cocktail glass model element by sharing the same nodes.

Governing Equations
Following Zienkiewicz and Bettss (1982), the governing equations for porous media saturated with pore water are given as a combination of equilibrium equation and mass balance equation of pore water as follows (u-p formulation).
where σ ij : stress tensor, k ij : coefficient of permeability, ε ij : strain tensor, ρ: mass density (saturated density), ρ f : mass density of pore water, g i : gravity acceleration vector, u i : displacement vector, n: porosity, K f : bulk modulus of pore water.

Boundary Conditions
The following boundary conditions are assigned.
where u i -Displacement specified on boundary Γ 1 . T i -Traction specified on boundary Γ 2 . p-Pore water pressure specified on Γ 3 . q-Pore water inflow specified on Γ 4 . _ w j -Velocity of pore water relative to soil skeleton (average over total sectional area).
The union of the boundaries Γ 1 and Γ 2 is equal to the total area Γ of the analysis domain V a . as shown in Fig. 25.2. The union of the boundaries Γ 3 and Γ 4 is equal to the total area Γ of the analysis domain V b .

Overview
The characterization tests, including cyclic undrained triaxial tests, of Ottawa F65 sand were conducted (El Ghoraiby et al. 2017. We performed the calibration of the model parameters for the constitutive model, that is, the cocktail glass model, as a part of the LEAP 2017 simulation exercise . We are concerned here with the calibration procedures and the model parameters obtained. Also, we will show a comparison between the results of the numerical element simulations and the results of the corresponding cyclic undrained triaxial tests.

Calibration Process (First Step)
We decided on the values of the parameters that can be determined without trial and error. Table 25.2 shows model parameter values determined through the first step. These parameter values were obtained through the following process: 1. Saturated mass density ρ sat was obtained from void ratio e 0 , specific gravity ρ s , and mass density of pore water ρ w . 2. Porosity n was calculated from void ratio e 0 . 3. Mean effective confining pressure p a corresponding to G ma , K L/Ua was set to 100 kPa. This is an arbitrary value. 4. Initial shear stiffness G ma for mean effective confining pressure p a was calculated using the following equation (Hardin and Richart Jr 1963).
Poisson's ratio ν was assumed to be 0.33. 6. Bulk modulus K L/Ua was calculated by the relationship between bulk modulus, initial shear modulus, and Poisson's ratio for isotropic elastic medium. 7. Using the following equation, internal friction angle ϕ f was calculated from the gradient M measured on the compression side of the effective stress path diagram.
8. Phase transformation angle ϕ p was assumed to be 28 degrees. 9. Upper bound for hysteretic damping factor h max was assumed to be 0.24.

Calibration Process (Second Step)
We decided on the values of the parameters that should be determined by trial and error. During the trial and error procedure, we kept the values determined through the first step. Table 25.3 shows the model parameter values determined through the According to the characterization test results, the axial strain time histories were unidirectional, but this tendency was not seen in the element simulations. Therefore, we did not select the liquefaction resistance curves determined with 2.5% single amplitude axial strain, shown in the report (El Ghoraiby et al. 2017, as targets of the calibration. Instead, we adjusted the parameters by fitting mainly axial strain double amplitude shown in the analytical time histories to those obtained from the triaxial tests.  Figure 25.5 shows sample comparisons of the effective stress path, the deviatoric stress-axial strain relation, the mean effective stress-axial strain relation and the time history of the excess pore water pressure ratio. Figure 25.6 shows comparisons of the time histories of mean effective stress and axial strain. In   Figure 25.7 shows details of the finite element model. We divided the whole area into four sub-areas so as to be able to assign the different material parameters to each area based on the soil density. However, we set the same material parameters to each area as discussed later. Figure 25.8 shows the finite element mesh and the boundary conditions. Each element shown in Fig. 25.8 is a cocktail glass model element. Every cocktail glass model element is accompanied by a pore water element (partially drained) by sharing the same nodal points. Total number of the cocktail glass model elements was 2686; hence, total number of the pore water elements was also 2686. Total number of the nodal points was 2800.

Model Geometry/Mesh and Boundary Conditions
The boundary conditions we applied are shown in Fig. 25.8.

Time integration method.
The coupled set of equations consisting of equilibrium equation and mass balance equation of pore water was integrated in time by using generalized time integration method (SSpj method, Zienkiewicz et al. (2005)). We set the time integration interval Δt to 0.001 sec so as to avoid the numerical instability in all cases. We also performed numerical analyses for 50 s, therefore, the number of time steps was 50,000 in each case.
where C is a Rayleigh damping matrix, M is a mass matrix. K is an initial tangential stiffness matrix at the beginning of dynamic analysis and is constant throughout the entire analysis time. We set parameters α, β as follows in all cases.
α ¼ 0:0, β ¼ 1:555 Â 10 À4 ð25:17Þ We performed eigenmode analysis under undrained condition in each centrifuge experimental case. The frequency of the first eigenmode was almost the same throughout the entire cases. They were approximately 10.24 Hz. Coefficient β in Eq. (25.16) was set to 1.555 Â 10 À4 so that the damping constant h of the first eigenmode becomes 0.005.

Material Properties and Constitutive Model Parameters
Model parameters determined in Phase I were used as they were, except for mass density, porosity and permeability coefficient. According to the following procedure, we calculated the dry density (ρ d ) of the ground for each case of the centrifuge model test.
1. We read the cone penetration resistance (q c ) at three or four depths (1.5 m, 2.0 m, 2.5 m, 3.0 m) from the graph representing the depth dependence of q c submitted by each institution conducting centrifuge model test. 2. The dry density (ρ d ) at each of these depths was read from the correlation diagram (Manzari 2017) of the cone penetration resistance (q c ) to the dry density (ρ d ). 3. The dry densities obtained at these depths were averaged. Table 25.4 shows mean dry density and coefficient of variation of dry densities for each centrifuge model experiment. As shown in Table 25.4, the coefficients of variation are mainly smaller than or equal to 1.0%; therefore, we decided to apply these mean values to the entire soil layer in each centrifuge model test.
In the Phase I, the parameter sets of the cocktail glass model were determined for three cases of dry densities; ρ d ¼ 1665.6, 1712.6, 1744.2 kg/m 3 . Each centrifuge test case had to be assigned with one of these parameter sets having the nearest dry density, because the liquefaction characteristics of soil with any other dry density except those of the three cases in Phase I was unknown.  As shown in Fig. 25.9, the dry density of each centrifuge test was closest to the dry density of the loosest case in Phase I. Therefore, in the type B simulation, we applied a parameter set of the loosest specimen, whose ρ d was 1665.6 kg/m 3 , to any case in the centrifuge model test. This implies that the liquefaction strength in all centrifuge model test cases becomes the same with each other, and then the results of the seismic response analyses in the type B simulations may become similar.
As shown in Table 25.4, saturated mass density (ρ sat ), porosity (n) and coefficient of permeability at T ¼ 20 C (kÃ) were determined for each centrifugal test case based on their mean value of mass density ρ d . Figures 25.10,25.11,25.12,25.13,25.14,and 25.15 show the contour plots obtained by initial self-weight analysis for the model of UCD1. Since the plot diagrams of the other cases are almost the same as the diagrams for UCD1, these drawings of the other cases are omitted. In these figures, compression of stress is represented by negative sign, and tension is represented by positive sign. Figures 25.10 and 25.11 show the distributions of vertical displacement and horizontal displacement, respectively. Figure 25.12 shows the distribution of pore  Initial state of the model for UCD1-shear stress τ xy (kPa) due to gravity water pressure. According to this figure, if the depth from the assumed water surface (see Fig. 25.8) is the same, the pore water pressure will be the same. The effective vertical stress distribution shown in Fig. 25

Results of the Dynamic Analysis
Seismic response analysis of 50 s was conducted for each case. The excitation was carried out in the first half for 20 s, and in the latter half for 30 s the drainage of the increased excess pore water pressure was tracked. Here, we show comparisons of time history between the experiment results and numerical simulation results.  Figure 25.19 shows the residual horizontal displacement distribution of the case "UCD3" obtained by numerical simulation. Figure 25.20 shows the comparisons of excess pore water pressure time histories at point P1. Figure 25.21 shows the excess pore water pressure distribution at time 15.1 s obtained by numerical simulation for the case "UCD3."  Figure 25.17 shows that the acceleration time histories by numerical simulations trace the results of experiments well, except spikes, which can be seen in the results of numerical simulations. These spikes seem to indicate that there is attenuation that is not taken into consideration in the analysis. For example, attenuation due to sliding friction between sand and rigid wall can be mentioned. In any case, if there is attenuation that is not considered in the analysis, the analysis result will have a larger displacement. Figure 25.18 shows that the horizontal displacement time histories by numerical simulations trace well the results of experiments in case of Ehime2, KAIST2, KyU3, and UCD3 except oscillations, which can be seen in the simulation results. The cause of these oscillations seems to be the same as the cause of the spikes in the acceleration time histories. If we take into account the attenuation, then the amplitude of the oscillations will decrease.
On the other hand, according to Fig. 25.9, the results of the numerical simulation may be underestimated for cases of "CU2", "KAIST2," "UCD3" and "ZJU2," and expected to be correct in other cases. Because, the dry densities of these four cases are below the dry density of the loosest sand in Phase I and the dry densities of the other cases are roughly equal to the dry density of the loosest sand in Phase I. If we take into account the unconsidered damping in the numerical analyses, then the residual displacements in Fig. 25.18 will decrease and the tendency of these residual displacements appears to become close to the above expectation based on Fig. 25.9.
According to the comparison of the time history of the excess pore water pressure shown in Fig. 25.20, the results of numerical simulation generally correlate with the test results except for the water pressure dissipation phase. However, in terms of maximum value of excess pore water pressure, the numerical simulation results in many cases exceed the test results.

Conclusion
As a part of the LEAP-UCD-2017 project, type B simulations were performed for centrifuge model tests using a two-dimensional effective stress analysis program "FLIP ROSE" incorporating a cocktail glass model (Iai et al. 2011). The centrifuge model tests were conducted as a part of the LEAP-UCD-2017 project to investigate the behavior of the liquefiable ground with gentle slope.
In Phase I of the numerical simulation, three parameter sets of the cocktail glass model were determined by trial and error so as to reproduce the results of the  undrained cyclic triaxial tests for the specimens of three different dry densities. In Phase II, models were made based on the data of the centrifuge model tests. For each centrifuge model test, we decided which parameter set of the cocktail glass model to apply by referring to the dry density of the ground. As a result, for all of the centrifuge model tests, the parameter set obtained based on the triaxial test result of the specimen with the loosest dry density was applied. It seems that it is difficult to adjust the dry density of the centrifuge model and in many of the models the density was lower than that of the loosest specimen.
Since it was a type B simulation, numerical simulations were performed without knowing the experiment results. Acceleration time histories, displacement time histories and excess pore water pressure time histories simulation results were generally consistent with the results of the centrifugal model tests. There was a tendency that the response values of the numerical simulation were slightly larger than the response values expected from the difference between the dry density of the ground of the centrifuge model and the dry density of the specimen on which the applied parameter set depends. These differences seem to be related to the attenuation due to sliding friction between the ground and the rigid wall, which is not taken into account in the numerical simulation.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.