Semi-Automated Procedure to Estimate Nonlinear Kinematic Hardening Model to Simulate the Nonlinear Dynamic Properties of Soil and Rock

: The strain-dependent nonlinear properties of ground materials, such as shear modulus degradation (G/G max ) and damping, are of signiﬁcant importance in seismic-related analyses. How-ever, the ABAQUS program lacks a comprehensive procedure to estimate parameters for a built-in model. In this study, a nonlinear kinematic hardening (NKH) model with three back-stress values was used, which allows better ﬁtting to the backbone curves compared to the simpliﬁed nonlinear kinematic hardening (SNKH) model previously proposed. Instead of modeling in ABAQUS, a semi-automated procedure was implemented in MATLAB, which can predict shear stress–shear strain hysteretic loops, to ﬁnd the ﬁtting parameters to the target G/G max and/or damping curves. The procedure was applied for three soil and two rock samples, and the results indicate a good match between model and target backbone curves, which proves the application of the procedure and the NKH model in simulating the nonlinear properties of ground materials.


Introduction
Shear modulus degradation (G/G max ) and damping (D) curves are the most important properties that need to be considered when simulating dynamic problems, such as seismic wave propagation in the ground, or dynamic soil-structure interactions (SSIs) [1][2][3][4]. Previous studies have indicated that good modeling of the nonlinear hysteretic properties of soil, including shear modulus degradation and damping, could enable wave propagation and soil-structure interaction to be accurately simulated [2,5].
Among famous commercial programs, such as ABAQUS, ANSYS, LS-DYNA, FLAC, and PLAXIS, ABAQUS has commonly been used for dynamic problems and SSI investigation, because of its capacity and convenience with numerous materials and elements for soil, structure, composite, and porous media. However, constitutive models with the implementation of G/G max and damping curves for soils and rocks in ABAQUS have still been quite limited. The simplified nonlinear kinematic hardening (SNKH) model was implemented in ABAQUS to simulate the strain-dependent properties of soil [5][6][7]. Although the procedure for estimating SNKH model parameters was proposed by Anastasopoulos and Gelagoti [5], with a single calibrated parameter based on the ultimate shear strength of material, the kinematic hardening parameters were varied to fit with G/G max obtained from experiments, without changing the shear strength at extremely large shear strains. Furthermore, the SNKH model did not accurately simulate the nonlinear dynamic properties of soil, due to the use of a simple nonlinear kinematic hardening rule. Several user-defined materials (e.g., UMAT and VUMAT) imply that more complicated hardening rules had been developed [4,[8][9][10]. However, the implementation of these UMAT materials A total of three soils were simulated in the study. Two soils, namely, a silica and a mixture of weathered soil and pebbles, were tested with the resonant column test; another popular soil, Toyoura sand, was also used in this study with the fundamental properties reported by Wu, Yamamoto [11] and Kokusho [12], as shown in Table 1. The samples prepared in laboratory RCT tests were made with a diameter of 50 mm and height of 105 mm. The silica sand was carefully poured into the cylindrical sampler at a relative density of 80%; the mixture was carefully tamped at an optimum water content, a wopt of 10.5%, defined from compaction test to obtain the highest density, as well as stiffness. The mixed soil included weathered soil and small pebbles with a mixed ratio by weight of 55% and 45%, respectively. Table 1 shows the strength properties of the two soils at similar conditions to the RCT test. The small-strain shear wave velocity, shear modulus degradation (G/Gmax), and damping (D) of the two soils were defined by performing resonant column tests (RCTs) at a confining pressure of 200 kPa, an average value of confining strain when performing RCTs. Figure 1a,b show G/Gmax and the damping curves of the two soils. Figure 1c shows the G/Gmax and D curves of the Toyoura sand, which were reported by Wichtmann and Triantafyllidis [13].

Torsional Shear Test Modeling
To obtain the G/G max and damping curves by simulation, a torsional shear (TS) test was modelled in ABAQUS. A cylindrical column was made with the same dimension as the TS test column in the laboratory, as shown in Figure 2a. The end of the column was fixed, and a torsional shear deformation was applied at the top-free end of the column as a torsional rotation with a sinusoidal form (Figure 2b). The torsional strain for a cylindrical column specimen could be defined based on the applied rotational angle (θ), radius (r), and height (l) of the cylindrical column, as follows [14,15]: Appl. Sci. 2021, 11, 8611 4 of 16 the experimental RCT. The peak amplitude of the torsional shear deformation applied at the top-free end was changed in ten levels to obtain the shear strain of the column, which was defined from the maximum rocking angle, in the range of 0.0001% to ~0.1%. The torsional shear stresses of elements located outside the column were captured, and the hysteretic loop of shear stress-shear strain could be obtained for each strain level. The G/Gmax and damping value could be easily defined from the measured hysteretic loops.

Simplified Nonlinear Kinematic Hardening (SNKH) Model
The nonlinear kinematic hardening model has frequently been used for soil in dynamic analyses using ABAQUS [5][6][7]16,17]. In previous studies, a simplified nonlinear kinematic hardening (SNKH) model was applied with the implied von-Mises criterion and a back stress (α), which control the translation of yield stress. The yield criterion is defined as: where σo is the initial yield stress at zero plastic strain. The increment of back stress is described as follows [5,18]: where ̅ is the increment of equivalent plastic strain. The three parameters required for the SNKH model are the initial yield stress (σo), initial kinematic hardening modulus (C1 is equal to Young's modulus E), and rate of decrement of kinematic hardening modulus (γ1). The model parameters could be estimated based on the assumption of equaling the stress at large plastic strain to the maximum yield stress (σy) of soil, which can be estimated from the soil strength parameters (i.e., friction angle ϕ, and cohesion, C) [5,7]. The σo value can be estimated as a fraction (λ) of the yield stress (σo = λ·σy), and the γ1 parameter could be calculated as: In this study, the equivalent radius r eq was taken to be 0.82 times greater than the radius of the cylindrical column [15]. The confining pressure was applied by surface pressure load and maintained during modelling to simulate a confining stress. The value of the confining pressure was chosen to be 200 kPa to produce a similar stress condition to the experimental RCT. The peak amplitude of the torsional shear deformation applied at the top-free end was changed in ten levels to obtain the shear strain of the column, which was defined from the maximum rocking angle, in the range of 0.0001% to~0.1%. The torsional shear stresses of elements located outside the column were captured, and the hysteretic loop of shear stress-shear strain could be obtained for each strain level. The G/G max and damping value could be easily defined from the measured hysteretic loops.

Simplified Nonlinear Kinematic Hardening (SNKH) Model
The nonlinear kinematic hardening model has frequently been used for soil in dynamic analyses using ABAQUS [5][6][7]16,17]. In previous studies, a simplified nonlinear kinematic hardening (SNKH) model was applied with the implied von-Mises criterion and a back stress (α), which control the translation of yield stress. The yield criterion is defined as: where σ o is the initial yield stress at zero plastic strain. The increment of back stress is described as follows [5,18]: .
where . ε pl is the increment of equivalent plastic strain. The three parameters required for the SNKH model are the initial yield stress (σ o ), initial kinematic hardening modulus (C 1 is equal to Young's modulus E), and rate of decrement of kinematic hardening modulus (γ 1 ). The model parameters could be estimated based on the assumption of equaling the stress at large plastic strain to the maximum yield stress (σ y ) of soil, which can be estimated from the soil strength parameters (i.e., friction angle φ, and cohesion, C) [5,7]. The σ o value can be estimated as a fraction (λ) of the yield stress (σ o = λ·σ y ), and the γ 1 parameter could be calculated as: By calibrating the λ parameter, which typically ranges from 0.1 to 0.3 for soil, the model parameters could be estimated to fit to the G/G max curve. The mean-average-error (MAE) was used to estimate the difference between the model and target curves, as: where A m and A t are the G/G max values obtained from the SNKH model and target curve, respectively. Figure 1a-c show the G/G max and D curves with SNKH model compared to the target curves with the model parameters indicated in Table 1. The λ parameter for the three soils varied (0.15 to 0.3). Note that the main point of the calibrating method was modifying the back-stress curve by changing the σ o and γ 1 parameters, without changing the ultimate shear strength, which is represented by summation of the σ o and C 1 /γ 1 values [5]. In this study, the λ parameter was chosen as the value that produced the lowest MAE, as shown in Figure 1d.
The degradation of shear modulus and damping properties simulated in ABAQUS by the SNKH model partially matched the target curves. However, there were several error points: (1) the threshold shear strain (ε SNKH th ), the strain level that separates material behavior from elastic to nonlinear plastic that can be estimated as shear strain at G/G max of 0.99 [19,20] in the model was much higher than those in experiments (ε RCT th ); (2) the G/G max curve at small strain range was not well simulated; and (3) the hysteretic damping in the SNKH model at low high strain level was much smaller than in experiment, as shown in Figure 1a-c. Similar errors could also be seen in the results reported by Tsinidis and Pitilakis [6] and Tsinidis and Pitilakis [7].

Nonlinear Kinematic Hardening (NKH) Model
In SNKH model, the G/G max and damping curves were not well modelled; therefore, in this study, the nonlinear kinematic hardening (NKH) model, which was basically developed from the SNKH model, was recommended. The NKH model has been a very common model used in mechanical engineering, and several researchers have proposed the procedure to estimate model parameters for metal materials based on experimental results [21][22][23][24]. Using a similar yield surface function and the von-Mises failure criterion to the SNKH model, the NKH model was upgraded to three kinematic hardening moduli (i.e., C 1 , C 2 , and C 3 ), and three rates of decrement of kinematic hardening modulus (i.e., γ 1 , γ 2 , and γ 3 , respectively). With the initial yielding stress, σ o , the model had seven input parameters, resulting in difficulties in finding appropriate parameters for ground materials. The back stress that included the three components, α 1 , α 2 , and α 3 with their evolution could be estimated from each pair of kinematic parameters (i.e., C 1 -γ 1 , C 2 -γ 2 , and C 3 -γ 3 , respectively) by Equation (3). Figure 3a shows the back-stress curves for silica sand in the SNKH and NKH models with the plastic shear strain using the equation proposed by Mahmoudi and Pezeshki-Najafabadi [21]. In the NKH model, the total back stress is the summation of three backstress components. It can be seen that with the increasing γ p eq or shear strain, each back-stress component including α 1 , α 1 , and α 3 , respectively, sequentially prevail in the total back-stress (α in NKH model), and the total back-stress curve can be divided into three segments by the limitation of each back stress component. Therefore, each pair of C i and γ i parameters controls the G/G max curve at each level of shear strain, as indicated in Figure 3b. A manual process was performed to estimate the seven input parameters for the NKH model. Firstly, the initial yield stress, σ o , was estimated from G and the threshold shear strain (ε RCT th ), which is defined from the G/G max curve. Compared to those in the SNKH model, the σ o value in the NKH model was slightly reduced, because the SNKH model produced a larger threshold strain. Due to the reduction in the σ o , three pairs of kinematic hardening parameters were firstly chosen to have slightly higher back stress at all strain range ( Figure 3a). For each trial, the G/G max curve was obtained, and each pair of kinematic hardening parameters was modified following the instruction indicated in Figure 3b. By variation in its slope, the target G/G max curves are divided into three parts, which are controlled by each pair of kinematic parameters, C 1 -γ 1 , C 2 -γ 2 , and C 3 -γ 3 , respectively. For each part of the G/G max curve with its pair of kinematic hardening parameters, the reduction in C and increment in γ resulted in the increase in slope of the G/G max curve, and vice versa, as shown in Figure 3b with blue and green arrows. This is because the reduction in C and/or increase in γ in each segment results in reduction in limitation of back stress; therefore, G reduces in that segment.
order from C1-γ1 to C2-γ2, and finally C3-γ3; and as a result, the G/Gmax curve could gradually be fitted to the experimental curve. Figure 4 illustrates the results of the nonlinear curve by using the NKH model with parameters for the three soils obtained after approximately 20 trials with a considerable time consumption. The NKH model perfectly simulated the G/Gmax curve at all ranges of shear strain level. In particular, the results for damping obtained from the NKH model were significantly improved compared to the SNKH model in Figure 1a,b, due to the reasonable threshold shear strain, and the back-stress curves with more kinematic hardening parameters. Anastasopoulos and Gelagoti [5] reported the requirements for the yield stress of material by controlling back-stress at large strain, which were implemented in the SNKH model. Figure 3a typically shows the backstress curve of silica sand in the NKH model, and at extremely large shear strain, the backstress value was slightly higher than that in the SNKH model, because of the small reduction in the initial yield stress, σo. Therefore, the requirement of the ultimate shear strength of soil could be allowed. Back stress,  The estimation process for kinematic hardening parameters should be performed in order from C 1 -γ 1 to C 2 -γ 2 , and finally C 3 -γ 3 ; and as a result, the G/G max curve could gradually be fitted to the experimental curve. Figure 4 illustrates the results of the nonlinear curve by using the NKH model with parameters for the three soils obtained after approximately 20 trials with a considerable time consumption. The NKH model perfectly simulated the G/G max curve at all ranges of shear strain level. In particular, the results for damping obtained from the NKH model were significantly improved compared to the SNKH model in Figure 1a,b, due to the reasonable threshold shear strain, and the backstress curves with more kinematic hardening parameters. Anastasopoulos and Gelagoti [5] reported the requirements for the yield stress of material by controlling back-stress at large strain, which were implemented in the SNKH model. Figure 3a typically shows the back-stress curve of silica sand in the NKH model, and at extremely large shear strain, the back-stress value was slightly higher than that in the SNKH model, because of the small reduction in the initial yield stress, σ o . Therefore, the requirement of the ultimate shear strength of soil could be allowed.   Figure 5 shows the algorithm of the program based on the ABAQUS calculation procedures [18]. For a given incremental shear strain Δεk, the von-Mises yield criterion was initially checked with the original predicted elastic stress (Δσ pr k) by the assumption of pure elastic  Figure 5 shows the algorithm of the program based on the ABAQUS calculation procedures [18]. For a given incremental shear strain ∆ε k , the von-Mises yield criterion was initially checked with the original predicted elastic stress (∆σ pr k ) by the assumption of pure elastic strain. Then, the iteration process was applied to determine the correct incremental plastic strain (∆ε p ), incremental back stress (∆α), and incremental total stress (∆σ) at a step.

MATLAB Algorithm Estimating Shear Stress-Shear Strain Loop
In the iteration process, new plastic strain was estimated based on the flow rule [18], and the back-stress values either for the positive (α p ) or negative (α n ) phase could be calculated as follows [21]: where α ip and α in are the maximum and minimum values of the back-stress component during positive and negative loading, respectively, and ∆ε p is the incremental plastic strain at the step. The total shear stress could be calculated from the incremental elastic strain (∆ε e ), which was calculated from the incremental strain and plastic strain at that step (∆ε e = ∆ε k − ∆ε p ), and the yield function, F, could be defined following Equation (2). The process of iteration finishes when F reaches the tolerated value [18].
where αip and αin are the maximum and minimum values of the back-stress component during positive and negative loading, respectively, and Δε p is the incremental plastic strain at the step. The total shear stress could be calculated from the incremental elastic strain (Δε e ), which was calculated from the incremental strain and plastic strain at that step (Δε e = Δεk − Δε p ), and the yield function, F, could be defined following Equation (2). The process of iteration finishes when F reaches the tolerated value [18]. The procedure was implemented in MATLAB, which can estimate shear stress-strain loops with a much faster speed compared to ABAQUS. The results obtained from this procedure were then validated by comparing the stress, hysteretic loops, G/Gmax, and damping curves obtained in the ABAQUS simulation. Figure 6a,b show a comparison of the stresses obtained in ABAQUS and MATLAB with 10 cycles of sinusoidal shear strains, and the identical results between hysteretic loops obtained in both ABAQUS and MATLAB were typical for the silica and Toyoura sands. The shear modulus (G) and damping (D) were then estimated and compared with the ABAQUS results for the silica and Toyoura sands, as shown in Figure 6c. Similar results were also obtained for the mixture soil. The good matching of the results in both the G/Gmax and D curves proves that the algorithm could predict the soil nonlinear properties in a wide range of shear strain of 0.0001% to 1%.  The procedure was implemented in MATLAB, which can estimate shear stress-strain loops with a much faster speed compared to ABAQUS. The results obtained from this procedure were then validated by comparing the stress, hysteretic loops, G/G max , and damping curves obtained in the ABAQUS simulation. Figure 6a,b show a comparison of the stresses obtained in ABAQUS and MATLAB with 10 cycles of sinusoidal shear strains, and the identical results between hysteretic loops obtained in both ABAQUS and MATLAB were typical for the silica and Toyoura sands. The shear modulus (G) and damping (D) were then estimated and compared with the ABAQUS results for the silica and Toyoura sands, as shown in Figure 6c. Similar results were also obtained for the mixture soil. The good matching of the results in both the G/G max and D curves proves that the algorithm could predict the soil nonlinear properties in a wide range of shear strain of 0.0001% to 1%.

Semi-Automated Procedure to Estimate the NKH Model Parameters
In the previous section, the NKH model parameters, including the initial yield stress σ o and three pairs of kinematic hardening parameters (i.e., C 1 -γ 1 , C 2 -γ 2 , and C 3 -γ 3 ), were found for the three soils by manually using the trial-and-error method, which consumed a lot of time and effort. Therefore, in this section, an algorithm is proposed with the main purpose of rapidly finding the NKH model parameters to match the target G/G max and damping curves obtained from experimental results based on the basic properties of soil, such as the Young's modulus (E) or shear modulus (G), density (ρ), and Poisson's ratio (ν). Figure 7 shows the flowchart of the semi-automated procedure estimating the NKH model parameters. The outputs of the process are the seven parameters of the NKH model previously mentioned from the input, including E or G, ρ, and ν, and the target G/G max and damping curves. The procedure could be divided into two parts: (1) finding the initial parameters from the G/G max curve; and (2) updating the parameters with consideration of both the G/G max and damping curves with a weighting factor (w).

Semi-Automated Procedure to Estimate the NKH Model Parameters
In the previous section, the NKH model parameters, including the initial yield stress σo and three pairs of kinematic hardening parameters (i.e., C1-γ1, C2-γ2, and C3-γ3), were found for the three soils by manually using the trial-and-error method, which consumed a lot of time and effort. Therefore, in this section, an algorithm is proposed with the main purpose of rapidly finding the NKH model parameters to match the target G/Gmax and also be calculated from the σ and σo values. Finally, by applying the equation of evolution of back stress with plastic strain (i.e., the positive phase of Equation (6)), the initial parameters (i.e., C1, γ1, C2, γ2, C3, and γ3) could be predicted. Figure 8 compares the G/Gmax and damping curves with the initial parameters with the target curves for the different materials. The results show slight differences between the G/Gmax curve obtained by the model with the target in the shear strain range (0.0005% to 0.01%), especially for the Toyoura sand in Figure 8c. Therefore, these initial parameters needed to be upgraded.  There are three steps to estimating the initial parameter from the G/G max curve. The first parameter, initial yield stress, σ o , was estimated from G and the threshold shear strain (ε RCT th ), as shown in Figure 7, following the method mentioned above. At a shear strain larger than ε RCT th , the G of the material was reduced compared to G max , and plastic strain and back stress are produced. In the SNKH and NKH models, the yield stress at a strain level is the summation of the initial yield stress and the back stress. Therefore, the backstress target curve with the plastic shear strain could be estimated for a target G/G max curve in the second step, as shown in Figure 7. The target G/G max curve was discretized to approximately 200 data points in the logarithm scale for the shear strain space. At the shear strain ε k with increment shear strain ∆ε k , the shear stress, σ, could be estimated by the shear modulus G calculated from the G/G max value. Plastic shear stress could be defined from the shear stress and the elastic shear stress, which is equal to σ/G max . Back stress could also be calculated from the σ and σ o values. Finally, by applying the equation of evolution of back stress with plastic strain (i.e., the positive phase of Equation (6)), the initial parameters (i.e., C 1 , γ 1 , C 2 , γ 2 , C 3 , and γ 3 ) could be predicted. Figure 8 compares the G/G max and damping curves with the initial parameters with the target curves for the different materials. The results show slight differences between the G/G max curve obtained by the model with the target in the shear strain range (0.0005% to 0.01%), especially for the Toyoura sand in Figure 8c. Therefore, these initial parameters needed to be upgraded.
In the second part of the procedure, the model parameters were modified to better fit the G/G max and/or damping curves. The G/G max and damping curves were divided into two or three segments by manually assigning the strain level at a slope-changing point, as shown in Figure 3b. Then, each pair of kinematic hardening parameters, including C 1 -γ 1 , C 2 -γ 2 , and C 3 -γ 3 , respectively, were modified with a reducing changing rate to fit each part of the G/G max and/or damping curves. The changing rate, which could be manually input into MATLAB, indicates the change in the percent of the parameter (i.e., 20%, 10%, and 5%), and the degree of change should depend on the MAE value. The weighting factor (w) was used to estimate the error considering both G/G max and damping with MAE G and MAE D , respectively, as follows:  In the second part of the procedure, the model parameters were modified to better fit the G/Gmax and/or damping curves. The G/Gmax and damping curves were divided into two or three segments by manually assigning the strain level at a slope-changing point, as shown in Figure 3b. Then, each pair of kinematic hardening parameters, including C1-γ1, C2-γ2, and C3-γ3, respectively, were modified with a reducing changing rate to fit each part of the G/Gmax and/or damping curves. The changing rate, which could be manually input into MATLAB, indicates the change in the percent of the parameter (i.e., 20%, 10%, and 5%), and the degree of change should depend on the MAE value. The weighting factor (w) was used to estimate the error considering both G/Gmax and damping with MAEG and MAED, respectively, as follows: Note that in fitting to the damping curve, the new target damping was considered by subtracting the damping (i.e., damping obtained from the RCT) to the initial damping value, which was defined as damping at the threshold strain. New back-stress parameters for each segment of G/G max or damping curves could be obtained with the lowest MAE value at that segment, following the method proposed by Anastasopoulos and Gelagoti [5]. The convergence of the procedure could be reached if the MAE value was stable, with slight modification in the rate of change. Figure 9a-c show the results of the G/G max and damping curves for the three soils with parameters obtained by the proposed procedure with w factors of 1, which means that only G/G max was considered during the modification process. Acceptable results were obtained in G/G max compared between the model and target curves, which show considerable improvement compared to the initial parameters in Figure 8a-c. By performing the algorithm in the MATLAB compiler, the results of the parameters were obtained in a very short time, compared to the manual trial-anderror method mentioned above. The results for the damping improved, even though the underrated damping was seen at small shear strain and damping was overestimated at large shear strains. Figure 9d compares the back-stress curves of Toyoura sand for three cases, including the SNKH model, the NKH model with parameters obtained from the semi-automated procedure (NKH_G), and from modifying the SNKH curves through trial-and-error (NKH_Trial). The back-stress value in large plastic strain value obtained by the semi-automated method is reasonable compared to those of the trial-and-error method and the SNKH model. The same results were also obtained for the silica sand and mixture soil. It could be concluded that the NKH model and its parameters obtained by the algorithm match the requirement for the shear strength of the material proposed by Anastasopoulos and Gelagoti [5].
Note that in fitting to the damping curve, the new target damping was considered by subtracting the damping (i.e., damping obtained from the RCT) to the initial damping value, which was defined as damping at the threshold strain. New back-stress parameters for each segment of G/Gmax or damping curves could be obtained with the lowest MAE value at that segment, following the method proposed by Anastasopoulos and Gelagoti [5]. The convergence of the procedure could be reached if the MAE value was stable, with slight modification in the rate of change. Figure 9a-c show the results of the G/Gmax and damping curves for the three soils with parameters obtained by the proposed procedure with w factors of 1, which means that only G/Gmax was considered during the modification process. Acceptable results were obtained in G/Gmax compared between the model and target curves, which show considerable improvement compared to the initial parameters in Figure 8a-c. By performing the algorithm in the MATLAB compiler, the results of the parameters were obtained in a very short time, compared to the manual trial-and-error method mentioned above. The results for the damping improved, even though the underrated damping was seen at small shear strain and damping was overestimated at large shear strains. Figure 9d compares the back-stress curves of Toyoura sand for three cases, including the SNKH model, the NKH model with parameters obtained from the semi-automated procedure (NKH_G), and from modifying the SNKH curves through trial-and-error (NKH_Trial). The back-stress value in large plastic strain value obtained by the semi-automated method is reasonable compared to those of the trial-and-error method and the SNKH model. The same results were also obtained for the silica sand and mixture soil. It could be concluded that the NKH model and its parameters obtained by the algorithm match the requirement for the shear strength of the material proposed by Anastasopoulos and Gelagoti [5].   Figure 10 shows the results of the shear modulus degradation and damping curves for the three soils with parameters obtained by the semi-automated procedure with a w factor of 0.95, which means that both the G/G max and damping were considered. The results for damping were much improved for all materials compared to those in Figure 9a-c with acceptable shear modulus degradation curves. The value of the weighting factor could be increased to improve the accuracy of the model in simulating the shear modulus of the material.

21, 11, x FOR PEER REVIEW
13 of 16 Figure 10 shows the results of the shear modulus degradation and damping curves for the three soils with parameters obtained by the semi-automated procedure with a w factor of 0.95, which means that both the G/Gmax and damping were considered. The results for damping were much improved for all materials compared to those in Figure 9a-c with acceptable shear modulus degradation curves. The value of the weighting factor could be increased to improve the accuracy of the model in simulating the shear modulus of the material.

Application of the NKH Model and the Semi-Automated Procedure for Rocks
Although the application of the nonlinear model for rock material during dynamic analyses has recently been seen frequently with the increase in the structure's size and weight, such as in nuclear power plants (NPPs), high-rise buildings, and gravity dams, a few nonlinear models did exist for rocks in ABAQUS. Drucker-Prager, an isotropic hardening model recommended for rock materials, requires yield stress and plastic strain data, which have to calibrated with complicated techniques [25]. Moreover, in rock engineering, the Young's modulus, shear wave velocity, Poisson's ratio, and rock quality are the common properties provided during the geotechnical survey, whereas the cohesion and friction angle are not very common. Nonlinear dynamic properties for intact or jointed rocks have been proposed in several reports [26][27][28][29]. EPRI [30] proposed referent G/G max and damping curves for rock material, depending on its depth. Therefore, by using the semi-automated procedure proposed above, the parameters for the NKH model can be found for rock based on common fundamental properties and referent dynamic curves. Table 2 shows the ground profile for Wolsung NPP located in Korea. There were 11 rock layers with shear wave velocities varying from about 549 m/s at the surface to around 2804 m/s at a depth of 100 m. The semi-automated procedure was applied to define the NKH model parameters for layers 3 and 5 at depths of 6 m and 17 m, respectively; based on these depths, the referent G/G max and damping curves are presented in EPRI [30]. Figure 11 shows the G/G max and damping curves for finding the NKH model for the two rocks, considering only the shear modulus (i.e., the weighting factor, w, is 1). Good results were obtained for G/G max for the two rocks, proving the application of the NKH model with the semi-automated procedure for both soil and rock materials.

Application of the NKH Model and the Semi-Automated Procedure for Rocks
Although the application of the nonlinear model for rock material during dynamic analyses has recently been seen frequently with the increase in the structure's size and weight, such as in nuclear power plants (NPPs), high-rise buildings, and gravity dams, a few nonlinear models did exist for rocks in ABAQUS. Drucker-Prager, an isotropic hardening model recommended for rock materials, requires yield stress and plastic strain data, which have to calibrated with complicated techniques [25]. Moreover, in rock engineering, the Young's modulus, shear wave velocity, Poisson's ratio, and rock quality are the common properties provided during the geotechnical survey, whereas the cohesion and friction angle are not very common. Nonlinear dynamic properties for intact or jointed rocks have been proposed in several reports [26][27][28][29]. EPRI [30] proposed referent G/Gmax and damping curves for rock material, depending on its depth. Therefore, by using the semiautomated procedure proposed above, the parameters for the NKH model can be found for rock based on common fundamental properties and referent dynamic curves. Table 2 shows the ground profile for Wolsung NPP located in Korea. There were 11 rock layers with shear wave velocities varying from about 549 m/s at the surface to around 2804 m/s at a depth of 100 m. The semi-automated procedure was applied to define the NKH model parameters for layers 3 and 5 at depths of 6 m and 17 m, respectively; based on these depths, the referent G/Gmax and damping curves are presented in EPRI [30]. Figure 11 shows the G/Gmax and damping curves for finding the NKH model for the two rocks, considering only the shear modulus (i.e., the weighting factor, w, is 1). Good results were obtained for G/Gmax for the two rocks, proving the application of the NKH model with the semi-automated procedure for both soil and rock materials.

Conclusions
The nonlinear kinematic hardening (NKH) model was used in this study to simulate the nonlinear strain-dependent dynamic properties of soils and rocks. Seven model parameters were estimated either by trial-and-error or by using a semi-automated procedure with the modeling of torsional shear test in either ABAQUS or MATLAB. The main conclusions are as follows:

•
The simplified nonlinear kinematic hardening (SNKH) model could simulate the strain-dependent nonlinear properties of soil with a simple procedure of calibrating model parameters using the yield strength of the material. However, limitations of the model were illustrated, such as the threshold shear strain, shear modulus, and damping at low strain level; • Parameters for the nonlinear kinematic hardening (NKH) model were estimated by using the proposed method of estimation based on the G/G max curve with initial values defined from the SNKH back-stress curve. The NKH model with seven parameters improved on the limitations of the SNKH model, and its simulated nonlinear shear modulus degradation and damping curves of soil were much better than the frequently used SNKH model; • The trial-and-error procedure was proposed to estimate NKH model parameters, with one drawback of considerable time consumed, which has resulted in the very limited implementation of the NKH model in geotechnical engineering; • A semi-automated procedure was proposed to estimate the NKH model parameters based on the basic properties of soil and the nonlinear dynamic curves. By implementation with the help of MATLAB, the parameters of the NKH model could easily be estimated in a very short time, and considering either the shear modulus degradation or the damping of material. The procedure was also applied for rock material, with good results being obtained; • In this study, the NKH model could accurately simulate the shear modulus degradation, but the damping was not well simulated, especially at the high shear strain level, which would be not appropriate for analyses of a strong earthquake. Moreover, a more automated procedure should be further developed in combination with ABAQUS for estimating NKH model parameters which would be more applicable for researchers and engineers. Additionally, a built-in model combined with yielding criteria for geotechnical material needs to be developed further for ABAQUS that has been available in other commercial programs.