Numerical Simulation Trial by Cocktail Glass Model in FLIP ROSE for LEAP-UCD-2017

LEAP-UCD-2017 is a blind prediction exercise for numerical modelers with relatively simple centrifuge tests. However, the application of numerical simulation should be discussed from both modeler ’ s and practitioner ’ s viewpoints. The authors applied one of the most common dynamic analysis programs in Japan (FLIP ROSE) for this exercise. The program, FLIP ROSE, has been continuously updated for 20 years with voluntary supports from practitioners. Using this scheme, the authors discuss several issues including the possible variation in parameter determination in practice and the effect of arti ﬁ cial Rayleigh damping on predictions. Although only a limited number of cases of centrifuge tests were analyzed, the application of a common analysis program by practitioners pointed out some issues to be addressed by numerical modelers.


Introduction
LEAP (Liquefaction Experiments and Analysis Projects) is an international joint research project to discuss centrifuge modeling and numerical modeling of liquefaction. The name LEAP was proposed by one of the authors at an early phase of the project. However, the authors did not participate in the experimental program since Hiroshima University does not have a centrifuge. In other word, the authors present a numerical simulation perspective. Furthermore, the authors' main concern is not developing numerical constitutive models, but how to use numerical simulation adequately in practice. In short, the authors belong to the practitioner's side.
LEAP will contribute to verification and validation of the dynamic analysis program of soils. This is quite important for practitioners. Currently, we often used dynamic analysis programs in practice. The reliability and possible degree of variation in the computed results are quite important issues for practitioners, and the participation of practitioners in simulation in LEAP will be beneficial to the project. Thus, the authors join the simulation exercise with one of the most common dynamic analysis program in Japan: FLIP ROSE (Finite element analysis program for LIquefaction Process, Response Of Soil-structure systems during Earthquakes).
FLIP ROSE is a FEM program with an advanced constitutive model of soil that considers liquefaction. Both the program and the constitutive model have been continuously updated for these 20 years. These updates have been supported by enthusiastic voluntary support by practitioners. This experience indicates the importance of the participation by practitioners. The trial of LEAP-UCD-2017 exercise explained here is within the same context. With the voluntary support of the practitioners related to FLIP ROSE (FLIP Consortium WG members), the authors discussed several issues including the possible variation in parameter determination in practice and the effect of artificial Rayleigh damping. These discussions may contribute to the future direction on exercise design.

Numerical Analysis Code and Model Features
FLIP ROSE (ver.7.3.0) with the cocktail glass model (Iai et al. 1992(Iai et al. , 2011 is used for element test simulation. The main features of the constitutive model are as follows: 1. The model in FLIP ROSE is based on multiple shear mechanisms. The contact forces between particles in an arbitrary direction are considered as the stress states of the soil element. 2. There are two types of constitutive models for liquefaction in FLIP ROSE. The cocktail glass model is the recent advanced model. In the cocktail glass model, dilatancy is given as the sum of contractive part and dilative part. The increment of contractive part of dilatancy ε c d is given by the strain increment in a multiple shear mechanism γ xy , and dilative part of dilatancy ε d d is given by the current shear strain γ xy in an arbitrary direction. The value of dilative part of dilatancy was determined so as to make no work under shear. As a result, the void ratio decreases from the initial value e 0 to a possible limiting value e min , but the current void ratio of the soil element depends on the current strain level. Thus the shape of possible points of void ratio is similar to the shape of a cocktail glass as shown in Fig. 31.1. One of the most unique features of the code is that more than 100 groups are using the program in practice in Japan. Before the implementation of the cocktail glass model, which is used here, another model (only applicable for fully drained or fully undrained conditions) was implemented in the FLIP ROSE program. This model was quite stable and showed good performance in the analysis of the structures during the 1995 Kobe earthquake. Therefore, many practitioners use it in seismic design of structures. These practitioners gradually started using the cocktail glass model in practice.

Difference Between the Method Used in Practice and the Method of This Exercise
In practice, the in-situ information should be referenced as much as possible. It is very important to have parameters with physical meaning. For example, the internal friction angle shall be determined by CD or CUB tests of in-situ samples. Shear modulus shall be determined by PS logging results. However, there are only re-constituted samples in this project. Thus our standard procedure in parameter setting cannot be used. The authors determined the parameters from the reported laboratory test results and/or standard values often used in practice in case of insufficient information.
On the other hand, there are different types of parameters, which have no apparent physical meaning. For example, the parameters controlling the dilatancy characteristics are quite model-dependent. Although these parameters may have a conceptual background, there is no explicit procedure to determine the value of the parameters. In practice, practitioners often conduct the trial-and-error fitting process to obtain the appropriate value for these parameters. Liquefaction resistance curves are often used as the target of trial-and-error to obtain the parameters for dilatancy. In this case, cyclic tests were simulated in the program to obtain the computed liquefaction resistance curves. (Iai et al. 1990) Fig. 31.1 Schematic illustration of void ratio variation in cocktail glass model (Iai et al. 2011) In this exercise, the authors conducted the analysis in a two-dimensional plane strain condition. This is because 2D analysis is preferred in most of the practical cases for its simplicity. Therefore, parameter setting is also conducted in the two-dimensional plane strain condition in the simulation of the cyclic test of the soil element. In other words, a simple shear test or torsional shear test is more appropriate as the element test in this project. However, there were results of cyclic tri-axial tests, and the authors simulated these tests in the 2D condition. This constituted a difficulty in the application of our standard parameter setting procedures in this project.

Concept of Parameters and the Detailed Parameter
Setting Procedure The followings are steps used in the project: 1. Reference confining pressure, P a Parameter values are given at a certain value of confining stress called reference confining pressure. The reference confining pressure, P a , was determined from experimental test conditions (¼100 kPa). In practice, basically, it should be the in situ confining pressure. 2. Initial shear modulus, G ma In practice, it should be given by the PS logging result. In the test, bender elements or the resonance column test may be good to obtain the exact value. Without these test results, initial shear modulus, G ma , was calculated using the following equation (Kokusyo 1980): Here, e 0 : initial void ratio, σ 0 c : effective confining pressure. The equation is suitable for round-shaped particles. If the particle is not round shaped, a different equation shall be used. Considering the small value of friction angle in this case, the above equation was used. 3. Initial bulk modulus, K La and K Ua Initial bulk modulus, K La , for the loading process, and K Ua for the unloading process were calculated using the following equation: Here, ν: Poisson's ratio of soil skeleton and is assumed as 0.33.

Parameters mG, mK
mG is the parameter to control the confining pressure dependency of the shear modulus, and mK is the parameter to control the confining pressure dependency of the bulk modulus. mG and mK were set 0.5 as the general values. 5. Bulk density, ρ Bulk density, ρ, was determined from the experimental test condition. 6. Porosity, n Porosity, n, was calculated using the following equation: 7. Maximum damping coefficient, h max Maximum damping coefficient, h max , was set 0.24 as the general value of sand.
In practice, cyclic test results should be used to obtain these dynamic characteristics of the soil. It relates to the hysteresis damping of soil. When the effect of dilatancy is neglected, i.e. fully drained condition, a hyperbolic relationship is assumed for the shear stress-shear strain relationship. The extended Masing rule is applied to obtain the stress-strain in unloading and reloading processes. In the extension of the Masing rule in the model, the limit value of hysteresis damping in quite large strain is adjusted to agree with the input value of the maximum damping coefficient. 8. Soil strength parameters: cohesion, c, and internal friction angle, ϕ f Cohesion, c, and internal friction angle, ϕ f , are determined from laboratory test results (stress path shown in Fig. 31.2 as an example). The average value of the test results was used. 9. Permeability, k Permeability is determined from experimental results (¼0.01 cm/s). However, the process of pore water pressure dissipation and consolidation of liquefied soil were skipped in the computation exercises. This is because the deformation of simple soil layers' model without any clayey layers often stop immediately after the cessation of shaking as shown in the LEAP2014 exercise (Tobita et al. 2015), and the settlement due to consolidation is not significant (no more than 3-5%) as indicated in the experimental chart (Ishihara and Yoshimine 1992) in many cases. Thus, from a practical point of view, the authors are afraid that the computation process after shaking stopped is time consuming and less important. The authors assume that the effect of the permeability in this computation exercise is not significant and skip the detailed consideration in this exercise. 10. Phase transformation angle, ϕ p Phase transformation angle, ϕ p , is determined from laboratory test results (stress paths) as shown in example in Fig. 31.2. 11. Liquefaction parameters: other parameters for dilatancy Liquefaction parameters are necessary for successful numerical simulations of cyclic shear tests. The simulations were done for direct shear tests in a 2D plane strain condition, and the liquefaction strength curves were the main targets to be simulated. Furthermore, not only the liquefaction strength curves, but also the characteristics of strain accumulation rates with the increase in number of loading cycles were carefully considered.
Since the target of the centrifuge tests is an inclined slope, soil behaviors under initial shear stress may have important effects on the result. In other words, Kα effect may be important. As shown in Fig. 31.3, the cocktail glass model can consider Kα effect with appropriate choice in the parameters, which depend on the relative densities. However, due to insufficient test results of the soil behaviors under initial shear, the detailed consideration of Kα effect was skipped in the parameter determination in this exercise.

Parameter Setting Results Considering the Variation in Densities
The parameters for Ottawa F65 in loose (e 0 ¼ 0.585), medium (e 0 ¼ 0.542), and dense conditions (e 0 ¼ 0.515) were determined. The values of parameters are shown in Table 31.1 and in Table 31.2. The comparison of target liquefaction resistance

Parameter Setting Results by Multiple Practitioners
The program FLIP ROSE has been continuously updated for the last 20 years with voluntary supports from practitioners. Using this scheme, the authors tried to discuss the possible variation in parameter determination in practice. Currently, the FLIP Consortium organizes a working group (WG) to discuss the appropriate procedure of parameter determination in practice. The WG members are practitioners with the experience of using FLIP ROSE in practice. The authors asked them to determine the possible parameters independently. The results from nine practitioners from four companies were obtained.
The values of determined parameters for Ottawa F95 sand under medium dense conditions are shown in Table 31.3. The differences in the determined values were caused by the following reasons: 1. The methods to determine the initial shear stiffness were different. In reality, there are more variations in the procedure to determine the shear stiffness, i.e., from the result of small shaking prior to the main shaking, or from the result of inflight cone penetration in the centrifuge tests. 2. Regarding the method to determine the internal friction angle ϕ f , some practitioners determine the friction angle from the observed stress path, and other practitioners apply an empirical formula. 3. In the model, the behavior of soil is controlled by several parameters. There are multiple combinations of parameter set to successfully simulate the target behavior of soil. In other words, there is no theoretical background for a unique combination of parameter set to be obtained. 4. The degree of agreement to the target is different. Not only for the overall agreement, but the stress level and/or strain level, which were the focus of the practitioner, may be different. For example, Fig. 31.5 shows the variation of computed liquefaction resistance curves defined by single amplitude of strain 2.5%. Generally, the computed results are similar but there still remains slight differences.

Target Test Case and Parameter Determination
A total of nine centrifuge tests were reported as the target of Type B simulation. There are some differences in test conditions, but the most important difference is the soil density. As mentioned above, the parameters for loose (e 0 ¼ 0.585), medium dense (e 0 ¼ 0.542), and dense (e 0 ¼ 0.515) conditions were determined. However, the actual void ratio is far larger than the loose condition in the parameter determination except for one case (KAIST-1). It means the parameter determination for eight cases shall be extrapolations from three reported cases, and the reliability of the parameter is quite low. Note the situation above is related to the significant variations of the reported values of the maximum and minimum void ratios for Ottawa F-65 sand. The reason for these variations is not clear, but the authors guess the difference in the sand lot may result in the difference of soil behavior. Anyway, this variation resulted in a significant uncertainty on the achieved relative density of the soil. Considering the target tests were done at George Washington University (GWU), the relative densities for loose (e 0 ¼ 0.585), medium dense (e 0 ¼ 0.542), and dense (e 0 ¼ 0.515) conditions using the reported values of the maximum and minimum void ratios from GWU are 62%, 80%, and 91%, respectively. Figure 31.6 show the schematics of parameter determination with interpolation and extrapolation. The left figure is for the parameter r εdc . The variation is quite linear as a function of void ratio, and both interpolation and extrapolation work well. However, the right figure is for the parameter q 2 . The variation of this parameter is quite non-linear, and the extrapolated value for large void ratio is negative; nevertheless the value of the parameter is theoretically positive. Therefore, the authors picked up only the KAIST-1 case for the Type B simulation, and the exact value of the parameters are determined by linear interpolation.

Analysis Conditions
An average of horizontal accelerations is used as the input motion. The original data showed some fluctuations on the time stamp; the authors estimate the fluctuation is just because of a small error in digitalization (cancellation of significant unit). Time step for all data is assumed as 0.005 s. Although there are small spikes on the recorded motion, no artificial manipulation was added to the input motion. Since the vertical motion is smaller than the horizontal motion, the authors ignored it due to time limitation to submit the predictions. The computation stopped when the input motion ends, and no computation on the consolidation process was conducted. This is because of time limitation and no interest of the authors in this issue. The pore water pressure dissipation process can be computed by FLIP ROSE/cocktail glass model (please refer the results of other teams in the LEAP exercises). As for the mesh geometry, the same mesh indicated as an example in the guideline is used (a total of 2470 nodes and 4480 elements including pore water elements and water elements above the ground surface). The bottom and side boundaries are fixed. The undrained boundary was given for pore water.

Artificial Rayleigh Damping
For numerical stability, the authors applied artificial Rayleigh damping of very small value. Rayleigh damping is given as α[M] + β[K], and we usually assume α ¼ 0 for FLIP application in practice.
The remaining Rayleigh damping parameter β for the FLIP program is carefully determined by parametric studies in 1D analysis. The idea of determination of Rayleigh damping is to avoid large values, which artificially reduce the deformation (Fig. 31.7).
Following the practical scheme of FLIP application, the parametric analysis for the left side column as shown in Fig. 31.8 were conducted. Dynamic analysis with a non-liquefaction condition (not considering dilatancy effect) with various values of Rayleigh damping parameter β was done.
With a larger value in damping parameter, a smaller value of the computed displacement in response of soil layers will be given. These kinds of reductions in response can be regarded as the result of unnecessary artificial damping. Therefore, the threshold value not to minimize the computed displacement is chosen as the Rayleigh damping parameter in practice. Figure 31.10 is the results of the case with KAIST-1 input motion. From this result, the authors chose β ¼ 0.0002 (α ¼ 0) as the Rayleigh damping parameter.   Figure 31.9 shows the acceleration and excess pore pressure response of the left array. There was no significant difference between the results in left, center, and right-side arrays. Spiky peaks were observed at all depth locations, and more spikes were observed at the points closer to the ground surface. Pore water pressure increased in the beginning; however, strong reduction of pore water pressure due to positive dilatancy was observed after 10 s. Slight reduction of pore water pressure was observed after 18 seconds at all points, and this may be related with dissipation. Many spiky peaks were observed in acceleration time histories. These spiky peaks may be related to instability of the analysis mainly induced by significant change of the soil status due to dilatancy. In many cases, these spikes can be reduced with artificial Rayleigh damping. However, using quite a small value of Rayleigh damping parameter (β ¼ 0.0002, α ¼ 0) in this case, the spiky peaks remained.

Deformation Pattern
A slight level of lateral spread of the ground is observed. Figure 31.10 shows the deformation after shaking (t ¼ 24.55 s). with the direction to make the ground flat: left side goes down, and right side goes up. However, as shown in Fig. 31.10, the ground was not flat at the end of shaking. This may be due to the positive dilatancy to restrain the deformation. Note the deformation in pore water pressure dissipation process is not included in this result.

Possible Variation on Estimated Deformation
Figures 31.11 shows a summary of the deformation at the center in all numerical cases. With increase in relative density, both horizontal displacement and settlement decrease. This trend is quite reasonable. With the cases of various parameters determined by many practitioners, the results are quite scattered. Although the parameters are determined for the target density (e 0 ¼ 0.558: between medium: e 0 ¼ 0.542 to loose: e 0 ¼ 0.585), the results are less than 1/10 of the medium sand result to 10 times of the loose sand results. It indicates that the magnitude of the deformation can be controlled by parameter setting, and comparison of the magnitude of deformation between the test results and computed results is not meaningful. In other words, the magnitude of deformation is not only the issue of constitutive modeling but also of parameter setting. A more comprehensive discussion on parameter determination procedure including laboratory testing method, Kα effect, etc., is necessary. Another problem is unstable computation with the parameter sets C-1 to D-3. This may be induced by insufficient Rayleigh damping in computation. As mentioned above, Rayleigh damping is determined by 1D analysis of the target soil layers, and the value of the damping should be re-determined when the soil profile has been changed. However, in this exercise, this process was skipped and only the soil parameters are modified in the parametric study for simplicity. In other words, with the variation of the parameter for soil layers, computations were done with inappropriate Rayleigh damping, and it induced variations in the computed results. Careful discussion on the damping effect in the computation should be done. Note, not only the artificial damping but also the numerical damping effect in the time integration scheme should be examined.

Conclusion
This study describes the LEAP-UCD-2017 exercise with one of the most common dynamic analysis programs in Japan (FLIP ROSE). Although only a limited number of cases of centrifuge tests were analyzed, the demonstration of a practical application of a common analysis program to this exercise clarified some issues that should be addressed by numerical modelers. First, there is a possible variation in parameter determination in practice. Second, the possible variation of parameters can cause a large variation in computed deformations. Third, the effect of artificial Rayleigh damping should be carefully examined.
More detailed discussions on the variation of the results and comparison with the centrifuge test results remain for future study.