Investigation of a fractional derivative creep model of clay and its numerical implementation

Computers and Geotechnics 119 (2020) 103387

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

Investigation of a fractional derivative creep model of clay and its numerical implementation

T

Xu-Bing Xu, Zhen-Dong Cui

⁎

State Key Laboratory for Geomechanics and Deep Underground Engineering, School of Mechanics and Civil Engineering, China University of Mining and Technology, Xuzhou, Jiangsu 221116, PR China

ARTICLE INFO

ABSTRACT

Keywords: Creep Triaxial test Fractional derivative Elastic viscoplastic Numerical simulation Secondary development

The fractional calculus has been successfully applied to characterize the rheological property of elastic viscoplastic materials. However, soils were seldom involved in fractional derivative constitutive models, and the topic of triaxial creep tests is rarely discussed through fractional calculus. This study presents a fractional derivative model that describes the time-dependent behavior of Shanghai marine clay. The fractional derivative Merchant (FDM) model is established by introducing fractional element into the Merchant model. The comparison of the calculated values of the FDM model and the Merchant model shows that the FDM model can better describe the instant elastic deformation, attenuation creep and steady-state creep stages of the clay. The FDM model was implemented into a finite-element code in ABAQUS, which is used to simulate the creep tests. The simulations demonstrate that the development model has good ability to predict the time-dependent behavior of clay.

1. Introduction Marine reclamation land is an important way to solve the shortage of land resources, which has been carried out by the Netherlands, Singapore, Mexico, Japan, and Dubai [1–3]. However, the settlement of many reclamation projects will continue to develop in the years or even decades after the completion of the project [4,5]. For example, the first phase of Kansai International Airport in Japan was completed in 1994. By December 2012, the average settlement of the first phase of the project reached 12.9 m, which is far more than its expected settlement [6]. The time-dependent creep behavior of soil is a major problem in predicting long-term settlement, especially for large-scale land reclamation projects [7–10]. Element models are used to describe creep characteristics of rock and soil, which are composed of Hooker elastomer, Newton viscous body, Saint-Venant plastic body and other elements in series or in parallel. The representative models of element models are Maxwell model [11], Kelvin model [12] and Merchant model [13]. Zhao and Zhou [12] considered the correction equation on thermal gradient and used the generalized Kelvin model to describe the axial creep deformation. Research suggests that the axial creep deformation behavior of the frozen saturated clay can be represented by generalized Kelvin model. The modified Merchant model was used to describe visco elasticplastic deformation of soil by Ye et al. [13]. The research indicates

⁎

that the modified Merchant model has highly flexible in simulating the regional land subsidence problem with very complicated deformation characteristics of sediments. The researches show that many elements are connected in series and parallel to describe the complex mechanical characteristics of soils [14–16]. However, there are many parameters in the model, and the parameters have no clear physical meaning. Fractional derivative operators can reflect historical dependence and population relativity, which is exactly in line with the research needs of non-Newtonian fluid mechanics, complex viscoelastic mechanics and porous media mechanics [17–19]. Gemant [20] first proposed a rheological model based on fractional derivatives. Since then, extensive research has been carried out in this field. Fractional calculus has achieved great success in describing the rheological properties of viscoplastic materials [21]. Koeller [18] replaced Newton's dashpot in the classical Kelvin model and Maxwell model with an element of fractional derivative. The research shows that the constitutive relationship of viscoelastic body based on fractional derivative was consistent with the results of Rabotnov's viscoelastic theory. Cai et al. [22] applied fractional derivatives to viscoelastic models, and obtained the fractional derivative Maxwell and Kelvin models. Compared with the integer derivative model, the fractional derivative model has the advantages of fewer parameters, simpler expression and more accurate calculation results. Further analysis reveals that the model can describe the creep behavior of viscoelastic materials. A nonlinear incremental

Corresponding author. E-mail address: [email protected] (Z.-D. Cui).

https://doi.org/10.1016/j.compgeo.2019.103387 Received 23 October 2019; Received in revised form 8 December 2019; Accepted 10 December 2019 0266-352X/ © 2019 Elsevier Ltd. All rights reserved.

Computers and Geotechnics 119 (2020) 103387

X.-B. Xu and Z.-D. Cui

Nomenclature 3

q0 q p ε

e vep

effective confining pressure initial deviatoric stress deviatoric stress increment effective mean stresse total vertical strain

E1, E2 η α Γ

model based on a simple exponential function was first suggested by Yin et al. [23]. Three simulation schemes for modeling undrained triaxial tests, drained triaxial tests were established. The enhanced version of the model showed good performances in simulating triaxial tests on Toyoura sand under various loading conditions. In order to avoid excessive settlement after construction, the fractional derivative model was used to describe the creep characteristics of Nansha clay by Zhu et al. [24]. The research shows that the fractional calculus can obtain a good fitting effect with fewer parameters. Therefore, fewer parameters in the fractional derivative model is beneficial to the compilation and development of the model in numerical simulation. It is difficult to solve the problem by analytical method because of the non-linearity stress-strain relationship and the complexity of load and boundary conditions. The results of numerical analysis are one of the important techniques for geotechnical engineers to judge the problem [25–27]. Researches have shown that ABAQUS, a finite element software, contains abundant material models, element models [28,29]. The constitutive models in ABAQUS can solve many geotechnical problems, such as the static and dynamic problems, especially the nonlinear problems [29]. However, the constitutive models in ABAQUS comply with the classical soil mechanics constitutive model theory [26]. The stress-strain relationship of clay during unloading and reloading is considered to be elastic and will not produce plastic strain, but plastic strain will occur during reloading in fact. Therefore, the constitutive models in ABAQUS cannot effectively simulate the elastic viscoplastic mechanical properties of soft clay. It is necessary to compile and develop the constitutive model of soft clay in ABAQUS. The application of the developed constitutive model reflects the combination of test and numerical simulation. Based on the measured data, Chai et al. [30] carried out the finite element analysis and predicted the mechanical performance of the embankment. By comparing the predictions with the field data, it is concluded that simple models such as modified Cam clay, with reliable values of model parameters, can produce reasonable results for the performance of embankments on clayey deposit improved through prefabricated vertical drains for a relatively short time period. In this study, the triaxial tests were carried out under different confining pressures, which obtained the creep curves of clay under the same deviating stress. The experimental results were used to fit the parameters of the fractional derivative Merchant (FDM) model. The results show that the FDM model can describe the whole creep process more accurately. Then, the FDM model was compiled and developed by using the User-defined Material (UMAT) of ABAQUS. Following numerical simulation of the creep test with the developed model, a comparison of simulated results and experimental results indicates that the developed creep model is workable and can provide reference for creep calculation in engineering practice.

elastic strain elastic viscoplastic strain compression modulus viscosity coefficient fractional order deviatoric stress increment gamma function

mean stress p′ decreases with the increase of excess pore water pressure. Augustesen et al. [32] considered that the creep test is to measure the relationship between strain and time under the effective stress constant conditions. Therefore, the undrained creep test is not the actual creep test. In this study, the drainage creep tests of Shanghai marine clay were carried out. In order to compare with drainage creep test, an undrained creep test was carried out under the same conditions. 2.1. Sample preparation and experimental procedures The soil samples are the undisturbed soft clay, which were taken from a planned site in Shanghai. The physical and mechanical parameters of the samples were obtained by the tests, which are shown in the Table 1. It can be seen that the soil samples have high water content and compressibility. The creep tests of samples were carried out by using the triaxial apparatus of Global Digital Systems (GDS), UK. The schematic diagram of the triaxial equipment and the stress state diagram of the soil sample in pressure chamber are presented in Fig. 2. The tests design values were calculated according to the physical and mechanical parameters of the samples and the depth of the soils, which are shown in the Table 2. The samples were saturated and consolidated under confining pressure 120 kPa and backpressure 100 kPa conditions, and then anisotropic consolidation was carried out. The criterion for completion of consolidation is that the excess pore water pressure does not dissipate and strain rate is less than five thousandths per hour. When the consolidation was completed, the confining pressure remains unchanged, and a load was applied to the samples. And then, the excess pore water pressure curves and the creep curves were obtained. 2.2. Test Results and Discussions Fig. 3 (a) presents the axial strain of the drainage tests of the different confining pressures under the same deviating stress. According to the analysis of tests data, the strain curves can be divided into three stages, the instant elastic deformation, attenuation creep and steady-

2. Test Conditions and Results Triaxial creep test remains the axial and radial stresses constant and measures the variation of vertical strain with time. Triaxial creep tests are divided into drainage tests and undrained tests [31], as shown in Fig. 1. The drainage test remains the effective mean stress p′ and the deviatoric stress q constant. The undrained test remains the average total stress P and the deviation stress q constant, while the effective

Fig. 1. Sketch of triaxial creep tests under drained and undrained conditions. 2

Computers and Geotechnics 119 (2020) 103387

X.-B. Xu and Z.-D. Cui

under the same confining pressure. The strain curve develops rapidly and tends to be stable gradually under drained condition. Compared with the drain test, the strain of the undrained test initially develops slowly, but the trend of attenuation with time is smaller. Finally, the undrained strain curve exceeds the drainage strain curve. Fig. 5(b) shows the variation of the excess pore water pressure with time under different drained conditions. Under the undrained test conditions, the excess pore water pressure increases rapidly and remains stable after reaching a certain value. The excess pore water pressure cannot dissipate in time during construction, which will cause excessive settlement of foundation after construction. Therefore, the drainage consolidation [33–35] is used to reinforce soft clay foundation in the construction process.

Table 1 Physical and mechanical parameters of the samples. Water content w (%)

Density ρ (g/cm3)

Specific gravity ds

Void ratio e0

Liquid limit w L (%)

Plastic limit w P (%)

Plastic index IP (%)

35.5

1.73

2.67

1.08

36.70

21.30

15.40

state creep. Because of the small deviating stress, there is no the accelerated creep stage of the samples under the confining pressures. The instantaneous strains of soil samples occur after the instantaneous deviating stress is applied. With the increase of the confining pressures, the instantaneous strains become smaller and the strain difference of the soil samples become smaller and smaller. When the confining pressures are 84 kPa and 168 kPa, the instantaneous strains are 0.1819% and 0.0366% respectively. When the confining pressure doubles, the instantaneous strain decreases by about 6 times, which indicates that the change of the soil strains with confining pressures are non-linear. When confining pressure is 84 kPa and 168 kPa, the strains are 1.8824% and 0.6867 at 100 h, respectively. The confining pressure doubled and the strain increased to 2.7 times. Fig. 3 (b) shows the strain rates of the drainage tests of different confining pressures under the same deviating stress. The smaller the confining pressure is, the greater the strain rate becomes. When the confining pressure is 84 kPa, the strain rate is the largest after the instantaneous load. The larger the confining pressure is, the shorter time required to reach steady creep. Fig. 4 shows the variations of the excess pore water pressure with time under different confining pressures. The axial strain and the excess pore water pressure contrastive analysis found that the excess pore water pressure does not dissipate completely when the creep reaches the steady-state creep. The larger the confining pressure is, the smaller the excess pore water pressure becomes, and the faster the excess pore water pressure dissipates. In theory, the soil will not deform and the pore water will bear all the load increment when the saturated soil is subjected to the instantaneous load. However, the excess pore water pressures of the different confining pressures are different under the same instantaneous loading, which shows that the soil did not reach 100% saturation before the test. The results indicate that the creep characteristics of clay is significant with confining pressure being a critical factor. Fig. 5(a) describes the strain curves of different test conditions

3. Fractional elastic viscoplastic model and parameters analysis 3.1. Fractional elastic viscoplastic model Based on the test results, the FDM model is used to describe the creep characteristics of clay. The FDM model is constructed by replacing a Newtonian dashpot with a fractional derivative element in the Merchant model as shown in Fig. 6. The total strain ε is

=

e

+

(1)

vep

where e and vep are strains of the elastic and the elastic viscoplastic bodies, respectively. When consolidation is completed, the confining pressure remains unchanged and the main stress increases. We have

= e

=

e

=

(2)

vep

(3)

E1

= E2

vep

+

vep

t

, 0

1

(4)

where equals q. E1 and E2 are the compression modulus; η is the viscosity coefficient. Laplace transforms were applied to Eqs. (1), (2), (3) and (4), respectively. The constitutive equation is obtained in the transform domain.

Fig. 2. The schematic diagram of the triaxial equipment. 3

Computers and Geotechnics 119 (2020) 103387

X.-B. Xu and Z.-D. Cui

Table 2 The schemes of samples creep tests. Test conditions

Effective confining pressure

Drain

84 126 168 210 84

Undrain

¯=

1 1 + E1 E2 + s

3

(kPa)

Initial deviatoric stress q0 (kPa)

Deviatoric stress increment q (kPa)

Effective mean stresse p (kPa)

36 54 72 90 36

20 20 20 20 20

103 151 199 247 103

¯

Therefore, there are many parameters in the complex model, and the parameters have no clear physical meaning. Based on the above analysis, the FDM model is used to study the characteristics of the soft clay in this study. Compared with the Merchant model, there is an additional parameter α in the FDM model. In this study, the parameter α was analyzed. Fig. 8 shows the stain curves of the drainage test of the different fractional orders α under the confining pressure 84 kPa conditions. With the increase of the fractional order α, the axial strain develops faster and the time to reach the final settlement becomes shorter. The FDM model can be regarded as a linear elastomer when α is 0 and as the Merchant model when α is 1. The axial strain of viscoelastic materials (0 < α < 1) increase slowly with time, which neither remains unchanged like a linear elastomer (α = 0), nor increases like the Merchant model (α = 1). The above analysis shows that the FDM model is more applicable and can reflect the non-linear strain process well.

(5)

where ¯ is the Laplace transform of σ; ¯ is the Laplace transform of ε.The inverse Laplace transform and the generalized Mittag-Leffler function [36] were applied to Eq. (5). we obtain vep

( 1) k ( k+ + 1)

= k= 0

t E2

k+1

, 0

1, k= 0, 1, 2···

( + 1) = where Γ is the gamma function; ( ) = 0 e t t 1dt . Substituting Eqs. (3) and (6) into Eq. (1) leads to =

E1

+ k= 0

( 1) k ( k+ + 1)

t E2

( )

(6) and

k+1

, 0

1, k= 0, 1, 2··· (7)

The creep constitutive equation of the Merchant model [13] is

=

E1

+

E2

1

e

E2

t

3.2. Determination and analysis of the FDM model parameters Eq. (7) was programmed in MATLAB. The tests data of the drainage tests were fitted and the parameters of the FDM model were obtained. The tests data of the drainage tests and the FDM model simulation results are shown in Fig. 9, and the parameters of the FDM model are shown in Table 4. The parameters of the correlation coefficients R2 are more than 0.999, which shows that the FDM model can fit the strain curves of clay very well. The creep characteristics of clay were studied by Zhu et al. [24], and the research has shown that the fractional derivative model can achieve a good fitting effect by using fewer parameters. In this study, the FDM model uses 4 parameters to describe the instant elastic deformation, attenuation creep and steady-state creep stages of the clay. Fig. 10 illustrates the distribution of parameters with the confining pressures under the drainage tests conditions. The linear functions are used to fit the parameters. The compression modulus E1 and E2 show a good linear distribution with the confining pressures in Fig. 10(a) and 10(b). But for the viscosity coefficient η, the linear function cannot reflect the variations of the viscosity coefficient η with the confining

(8)

The creep constitutive equation of the Burgers model [14] is

=

E1

+

t+

E2

1

e

E2

t

(9)

In order to verify the practicability of the FDM model, the typical models, the Merchant model and the Burgers model, were used to compare with the fitness of the FDM model under the same confining pressure conditions. The results are shown in Table 3. Compared with the Merchant model, the Burgers model has a Newton dashpot. The correlation coefficient R2 of the Merchant model and the Burgers model are 0.9718 and 0.9959, which shows the fitness of the Burgers model is better than the Merchant model. But Compared with the typical models, the FDM model can get higher accuracy and better forecast effect. The same conclusion can be drawn from Fig. 7. Many elements are connected in series and parallel in many literatures [14–16] because the classical elements models cannot well reflect the characteristics of clay.

Fig. 3. Variation of (a) axial strain and time, (b) axial strain rate and time. 4

Computers and Geotechnics 119 (2020) 103387

X.-B. Xu and Z.-D. Cui

Fig. 4. Variations of excess pore water pressure and time (a)

3

pressures. The research suggests that the viscous coefficient reflects the thickness of the bound water film. When the soil is subjected to larger external loads, the thickness of the bound water film decreases. The stronger the bonding strength of water, the greater the resistance between soils. The stronger the viscosity of the bound water, the larger the viscosity coefficient, and then the final settlement time of the soil is enlarged. The results were put forward by Zhang et al. [37]. In Fig. 10(d), the results show the slope of the fitting function of the fractional order α is 0.00004, which indicates that fractional order α does little change with the confining pressures. Soft structured clays usually exhibit complex behaviors, which can lead to difficulties in the determination of parameters and high testing costs [38]. Besides, due to

= 84 kPa , (b)

3

= 126 kPa , (c)

3

= 168 kPa , (d)

3

= 210 kPa .

limited experimental resources, it is impossible to test all samples under different confining pressures conditions. The linear analysis of the model parameters can effectively predict model parameters. Fig. 11 shows the fitting results of the strain under different drainage tests conditions. It can be observed that both drainage and undrained tests can be well fitted by the FDM model, which shows that the applicability of the FDM model is very wide. The parameters of the FDM model under different drainage conditions are shown in Table 5. The compressive moduli E1 and E2 of the undrained test are larger than that of the drained test, which indicates that the excess pore water pressure is difficult to dissipate and soil is difficult to be compressed under undrained test conditions. However, the viscous coefficient η of

Fig. 5. Variations of (a) axial strain and time, (b) excess pore water pressure and time.

5

Computers and Geotechnics 119 (2020) 103387

X.-B. Xu and Z.-D. Cui

Fig. 6. Fractional derivative Merchant model.

the undrained test is smaller than that of the drained test, because the interaction force between soil particles decreases under undrained test conditions. Though, the fractional order α of the undrained condition is less than that of the drainage condition, the strain curve of the undrained test is finally greater than that under drainage test conditions. Further analysis reveals that the strain rate is not only affected by the fractional order α, but also affected by both the compressive moduli E1 and E2 and the viscous coefficient η.

Fig. 7. Comparison between the FDM model and the typical models predicted the test curve.

4. The realization of the FDM model in ABAQUS The finite element software ABAQUS has a wide range of applications and powerful functions, which can solve various problems from simple linear to complex non-linear. ABAQUS not only provides a standard finite element analysis program, but also has good open performance [25–27]. Users can use the secondary development platform to develop user subprogram interface and application program interface, and then the function of the main program is expanded, which has been widely used in practice [26–29]. The UMAT subroutine of the FDM model is written in Fortran language using the secondary development platform in this study. According to the requirements of the secondary development of ABAQUS, the developer needs to define Jacobian matrix of the integral point of the element material. The stiffness coefficient matrix of the constitutive relation of the material is

[D] =

Fig. 8. Axial stain curves of the drainage test of the different fractional orders a under the confining pressure 84 kPa.

(10)

where σ is the Cauchy stress tensor, i.e. the true stress; Δσ is stress increment; ε is strain tensor; Δ ε is strain increment. According to ABAQUS non-linear incremental loading technology and balanced iteration method, ABAQUS will call the user subroutine UMAT at the beginning of each incremental loading step, and calculate the stiffness coefficient matrix [D] by user program code. Then, ABAQUS forms the stiffness matrix [K]. The displacement increment ΔU is solved by the stiffness matrix [K] and the current load increment ΔR. The UMAT subroutine of the FDM model is shown in Appendix A. The stiffness matrix [K] of the FDM model is solved in Appendix A.

In order to verify the correctness of the UMAT subroutine, the drainage creep tests were simulated by ABAQUS. A cylindrical finite element model with a diameter of 50 mm and a height of 100 mm was established. The model was divided into 3 200 elements and the degrees of freedom in the x, y and z directions at the bottom of the model were constrained. The meshing of the finite element model of the sample is shown in Fig. 12(a). The calculation and analysis process were divided into two steps. The first step simulates the process of anisotropic consolidation, and the second step simulates the creep process in which the confining pressure remains constant and the

Table 3 The parameters of the models under the drainage test. Model

Effective confining pressure

FDM Merchant Burgers

84 84 84

3

(kPa)

E1 (kPa)

E2 (kPa)

109.95 109.95 109.95

59.14 12.45 16.13

(kPa·h) 10.09 149.06 103.27

6

(kPa·h) – – 3984.38

α

Correlation coefficient R2

0.6312 – –

0.9992 0.9718 0.9959

Computers and Geotechnics 119 (2020) 103387

X.-B. Xu and Z.-D. Cui

Fig. 11. Fitting results of the FDM model under different test conditions.

Fig. 9. Comparison between the tests data of the drainage tests and the FDM model predicted creep curves.

principal stress increases by 20 kPa. The model parameters are shown in Table 4. Fig. 12(b) shows the displacement picture when the confining pressure is 126 kPa. Fig. 13 illustrates the ABAQUS calculation results are basically consistent with the drainage tests data under different confining pressures. The UMAT subroutine of the FDM model can be effectively used to describe the instant elastic deformation, attenuation creep and steady-state creep stages of clay. It not only proves that the FDM model can describe the creep characteristics of clay well, but also proves the correctness of the UMAT subroutine of the FDM model in the secondary

Table 4 The simulated parameters of the FDM model. Effective confining pressure 3 (kPa)

E1 (kPa)

E2 (kPa)

η (kPa·h)

α

Correlation coefficient R2

84 126 168 210

109.95 371.75 546.45 706.71

59.14 142.52 211.87 284.18

10.09 16.86 23.89 23.73

0.6312 0.6234 0.6283 0.6349

0.9992 0.9995 0.9999 0.9999

Fig. 10. Variation of the FDM model parameters with confining pressure (a) E1, (b) E2, (3) η, (4) α. 7

Computers and Geotechnics 119 (2020) 103387

X.-B. Xu and Z.-D. Cui

Table 5 The parameters of the FDM model under different test conditions. Test conditions

Effective confining pressure

Drain Undrain

84 84

3

(kPa)

E1 (kPa)

E2 (kPa)

η (kPa·h)

α

Correlation coefficient R2

109.9505 140.1542

59.1422 62.6256

10.0914 4.4223

0.6312 0.4477

0.9992 0.9987

development of ABAQUS. The compilation and development of the FDM model is an effective combination of the laboratory tests and numerical simulation, which opens up a new field for the application of soils model. 5. Example and discussions From the above analysis, the FDM model can well describe the mechanical characteristics of clay. In this study, a foundation model of the Kansai International Airport was established, and the clay layer of consolidation settlement analysis was carried out with the UMAT subroutine of the FDM model and the sand layer of consolidation settlement analysis was carried out with the Mohr-colomb model. Fig. 14 [39] and Fig. 15 [4] describe the distribution of the soil layers on the seabed foundation. The parameters of the model were calculated according to the physical and mechanical parameters of the foundation, which were listed in reference [4]. The foundation model thickness is 100 m, and the length and width of the model are 700 m and 500 m as shown in Fig. 16. Before loading, the soil has been consolidated under the self-weight and the surface uniform load of 10 kPa. Then, according to the monitoring data, the load was simplified as a multi-segment linear load and applied to the surface of soil foundation. The monitoring data and the simplified results are presented in Fig. 17(a). In this study, the vertical compression of subseabed foundation was calculated by ABAQUS, and the results are shown in Fig. 16. In the subseabed foundation, the Ma13 was calculated by the UMAT subroutine of the FDM model. The thickness of Ma13 soil layer is 20.7 m. The vertical compression of Ma13 in numerical simulation was

Fig. 13. Comparison between the drainage tests data and the ABAQUS calculation results.

extracted and compared with monitoring data [4]. Fig. 17(b) illustrates the comparison between the monitoring data and the numerical calculation results. The results show that the FDM model can well reflect the compression characteristics of soft clay and predict the settlement of foundation. Fig. 12. Pictures of (a) Mesh diagram of computational model, (b) displacement under 3 = 126 kPa conditions.

8

Computers and Geotechnics 119 (2020) 103387

X.-B. Xu and Z.-D. Cui

ground settlement. The application of this model in the analysis of land subsidence is very meaningful. 6. Conclusions A series of tests for the Shanghai marine clay were conducted to investigate the compressibility behavior, creep characteristics after the instantaneous loading stage. The correlation of creep behaviors and the confining pressure of clay are interpreted. And then, the UMAT subroutine of the FDM model was used to simulate the creep tests, which enabled the combination of the laboratory test and numerical simulation effective. All these help us understand the long-term non-linear creep characteristics of the Shanghai marine clay. Based on the work presented and discussed, the following conclusions can be drawn. (a) The creep characteristics of clay is significant with confining pressure being a critical factor. With the increase of the confining pressure, the strain difference of the soil samples become smaller and smaller under the drainage conditions. The smaller the confining pressure is, the greater the strain rate becomes. The further study shows that the larger the confining pressure is, the smaller the excess pore water pressure changes, and the faster the excess pore water pressure dissipates. (b) Compared with the drain test, the strain curve of the undrained test initially develops slowly under the confining pressure 84 kPa conditions, but the trend of attenuation with time is smaller. And the excess pore water pressure increases rapidly and remains stable after reaching a certain value. (c) The FDM model simulates the instant elastic deformation, attenuation creep and steady-state creep stages of the clay with 4 parameters under the drainage conditions. The compression modulus E1 and E2 show a good linear distribution with the confining pressures. The linear function cannot respond well to the variations of the viscosity coefficient with confining pressure. The fractional order α does little change with the confining pressures. (d) The drainage and undrained tests data were fitted by the FDM model very well, which indicates that the applicability of the FDM model is very wide. The strain rate is not only affected by the fractional order α, but also affected by both the compressive moduli E1 and E2 and the viscous coefficient η. (e) The UMAT subroutine of the FDM model can well reflect the compression characteristics of soft clay and predict the settlement of foundation, which is effective combination of the laboratory tests and the numerical simulation.

Fig. 14. Subseabed Foundation Model.

Fig. 15. Subseabed profile at the Kansai International Airport site.

Declaration of Competing Interest

The further analysis indicates that the UMAT subroutine of the FDM model can be applied to analyze the site soil successfully. This study is a combination of the laboratory tests and the numerical simulation. If more soil samples can be obtained for experiment and numerical simulation, the calculation results can play a predictive role in the

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 16. Nephogram of the vertical compression.

9

Computers and Geotechnics 119 (2020) 103387

X.-B. Xu and Z.-D. Cui

Fig. 17. Load at seabed used in analysis together with observed and computed vertical compression of Ma13 (a) load at seabed, (b) vertical compression.

Acknowledgements

and Development Program (Grant No. 2017YFC1500702) and by the Fundamental Research Funds for the Central Universities (2015XKMS016).

This work in this paper was funded by the National Key Research Appendix A The major variables of the UMAT subroutine of the FDM model are

1. DDSDDT (NTENS): Jacobian matrix of the constitutive relations. 2. STATEV (NSTATV): An array is used to store state variables. The value is passed from the main program to the subroutine at the beginning of the incremental step, or updated and passed to the UMAT subroutine in the subroutine USDFLD or UEXPAN. At the end of the incremental step, the value of the STATEV must be returned to the main program. 3. NDI: Number of normal stress components. 4. NTENS: The size of the stress and strain component array (NTENS = NDI + NSHR). 5. NSTATV: Number of state variables, and it is defined by command * Depvar in ABAQUS. 6. PROPS (NPROPS): User-defined material parameters, that is, model parameters of user-defined material constitutive relationship. 7. NPROPS: Number of parameters of user-defined constitutive model. SUBROUTINE UMAT (STRESS, STATEV, DDSDDE, SSE, SPD, SCD, RPL, DDSDDT, 1 DRPLDE, DRPLDT, STRAN, DSTRAN, TIME, DTIME, TEMP, DTEMP, PREDEF, DPRED, 2 CMNAME, NDI, NSHR, NTENS, NSTATV, PROPS, NPROPS, COORDS, DROT, 3 PNEWDT, CELENT, DFGRD0, DFGRD1, NOEL, NPT, LAYER, KSPT, KSTEP, KINC) INCLUDE 'ABA_PARAM.INC' CHARACTER*80 CMNAME REAL*8 STRESS(NTENS), STATEV(NSTATV), DDSDDE (NTENS, NTENS), 1 DDSDDT(NTENS), DRPLDE(NTENS), STRAN(NTENS), DSTRAN(NTENS), 2 PREDEF (1), DPRED (1), PROPS(NPROPS), COORDS (3), DROT (3,3), 3 DFGRD0(3,3), DFGRD1(3,3), TIME (2) PARAMETER (zero = 0. d0, one = 1. d0, two = 2. d0, three = 3. d0, third = 1. d0 / 3. d0, half = 0.5d0) REAL*8 EMOD, ENU, EBULK3, EG2, ELAM REAL*8 GAMMA (2 0 1), DEEC (6), STRESS1(3) C Reading material parameters EMOD = PROPS (1) EMOD2 = PROPS (2) ETA = PROPS (3) ALAFA = PROPS (4) ENU = 0.4 EBULK3 = EMOD/(ONE-TWO*ENU) EG2 = EMOD/(ONE + ENU) EG = EG2/TWO EG3 = THREE*EG ELAM=(EBULK3-EG2)/THREE C Computational stiffness matrix DO K1 = 1, NDI DO K2 = 1, NDI DDSDDE (K2, K1) = ELAM END DO DDSDDE (K1, K1) = EG2 + ELAM END DO DO K1 = NDI + 1, NTENS DDSDDE (K1, K1) = EG END DO C Extracting gamma functions from TXT documents

10

Computers and Geotechnics 119 (2020) 103387

X.-B. Xu and Z.-D. Cui

C

C

C

GAMMA (:) = 0.0 open (unit = 151, file='D:\ gamma.txt') do I = 1,201 read (151, *, IOSTAT = istat), GAMMA(I) if ((istat /= 0)) then exit end if end do close (151) STATEV (11) = GAMMA (1) Extraction of Creep Stress State IF (KSTEP.EQ.2) THEN IF (TIME (1). LT. DTIME) THEN STRESS1(1) = STRESS (1) STRESS1(2) = STRESS (2) STRESS1(3) = STRESS (3) ENDIF STATEV (12) = STRESS1(3) STRESS1(1) = 0.0 STRESS1(2) = 0.0 STRESS1(3) = −0.02 Calculating creep strain increment A2 = 0.0 K2 = K3-1 IF (TIME (1). EQ.0) THEN DEEC (:) = 0.0 ELSEIF (TIME (1). EQ. DTIME) THEN DEEC (:) = 0.0 DO K1 = 1,3 DO K3 = 1,201 K2 = K3-1 A2 = A2 + STRESS1(K1)/ETA*(−1) **(K2)/GAMMA(K3) * * (ETA*TIME(1)**ALAFA/EMOD2)**(K2 + 1) END DO DEEC(K1) = A2 END DO ELSEIF (TIME (1). GT. DTIME) THEN DEEC (:) = 0.0 DO K1 = 1,3 A1 = 0.0 DO K3 = 1,201 K2 = K3-1 A1 = A1 + STRESS1(K1)/ETA*(−1) **(K2)/GAMMA(K3) * * (ETA*TIME(1)**ALAFA/EMOD2)**(K2 + 1) A2 = A2 + STRESS1(K1)/ETA*(−1) **(K2)/GAMMA(K3) * * (ETA*TIME(1)**ALAFA/EMOD2)**(K2 + 1)/ * TIME (1) *(ALAFA*(K2 + 1)) *Dtime ENDDO DEEC(K1) = A2 STATEV(6 + K1) = A1 END DO END IF STATEV (1) = DEEC (1) STATEV (2) = DEEC (2) STATEV (3) = DEEC (3) STATEV (4) = STATEV (4) + DEEC (1) STATEV (5) = STATEV (5) + DEEC (2) STATEV (6) = STATEV (6) + DEEC (3) DEEC (4) = 0.0 DEEC (5) = 0.0 DEEC (6) = 0.0 ELSE DEEC (:) = 0.0 END IF Calculating Cauchy stress DO K1 = 1, NTENS DO K2 = 1, NTENS STRESS (K2) = STRESS (K2) + DDSDDE (K2, K1) *(DSTRAN(K1)-DEEC(K1)) END DO END DO RETURN END

[2] Hoeksema RJ. Three stages in the history of land reclamation in the Netherlands. Irrig Drain 2007;56(S1):S113–26. [3] Wang W, Liu H, Li Y, Su J. Development and management of land reclamation in China. Ocean Coast Manage 2014;102:415–25. [4] Mesri G, Funk JR. Settlement of the Kansai International Airport Islands. J Geotech

References [1] Glaser R, Walsh PHPD. Land Reclamation in Singapore, Hong Kong and Macau. GeoJournal 1991;24(4):365–73.

11

Computers and Geotechnics 119 (2020) 103387

X.-B. Xu and Z.-D. Cui Geoenviron Eng 2015;141(2). [5] Furudoi T. The second phase construction of Kansai International Airport considering the large and long-term settlement of the clay deposits. Soils Found 2010;50(6):805-816. [6] Puzrin AM, Alonso EE, Pinyol NM. Unexpected excessive settlements: Kansai International Airport, Japan. Geomechanics of Failures. Dordrecht: Springer Netherlands; 2010. p. 23–43. [7] Graham J, Crooks JHA, Bell AL. Time effects on the stress-strain behaviour of natural soft clays; 1983;33(3):327–40. [8] Leroueil S, Kabbaj M, Tavenas F, Bouchard R. Stress–strain–strain rate relation for the compressibility of sensitive natural clays. 1985;35(2):159–80. [9] Yin JH, Graham J. Viscous–elastic–plastic modelling of one-dimensional time-dependent behaviour of clays. Can Geotech J 1989;26(2):199–209. [10] Zhu QY, Yin ZY, Hicher PY, Shen SL. Nonlinearity of one-dimensional creep characteristics of soft clays. Acta Geotech 2016;11(4):887–900. [11] Ter-Martirosyan ZG, Ter-Martirosyan AZ. Rheological properties of soil subject to shear. Soil Mech Found Eng 2013;49(6):219–26. [12] Zhao XD, Zhou GQ. Experimental study on the creep behavior of frozen clay with thermal gradient. Cold Reg Sci Technol 2013;86:127–32. [13] Ye SJ, Xue YQ, Wu JH, Li QF. Modeling visco-elastic–plastic deformation of soil with modified Merchant model. Environ Earth Sci 2012;66(5):1497–504. [14] Cong L, Hu X. Triaxial rheological property of sandstone under low confining pressure. Eng Geol 2017;231:45–55. [15] Zhao Y, Wang Y, Wang W, Wan W, Tang J. Modeling of non-linear rheological behavior of hard rock using triaxial rheological experiment. Int J of Rock Mech Min Sci 2017;93:66–75. [16] Tang LS, Zhao ZL, Chen Hk WuYP, Zeng YC. Dynamic stress accumulation model of granite residual soil under cyclic loading based on small-size creep tests. J Central South Univ 2019;26(3):728–42. [17] Bagley RL, Torvik PJ. A theoretical basis for the application of fractional calculus to viscoelasticity. J Rheol 1983;27(3):201–10. [18] Koeller RC. Applications of fractional calculus to the theory of viscoelasticity. J Appl Mech 1984;51(2):299–307. [19] Welch SJ, Rorrer RAL, Duren RG. Application of time-based fractional calculus methods to viscoelastic creep and stress relaxation of materials. Mech TimeDependent Mater 1999;3(3):279–303. [20] Gemant AXLV. On fractional differentials. London Edinburgh Dublin Philosoph Mag J Sci 1938;25(168):540–9. [21] Mainardi F. An historical perspective on fractional calculus in linear viscoelasticity. Fract Calculus Appl Anal 2012;15(4):712–7. [22] Cai W, Chen W, Xu W. Characterizing the creep of viscoelastic materials by fractal derivative models. Int J Non Linear Mech 2016;87:58–63. [23] Yin ZY, Wu ZX, Hicher PY. Modeling monotonic and cyclic behavior of granular

materials by exponential constitutive function. 2018;144(4):04018014. [24] Zhu HH, Zhang CC, Mei GX, Shi B, Gao L. Prediction of one-dimensional compression behavior of Nansha clay using fractional derivatives. Mar Georesour Geotechnol 2017;35(5):688–97. [25] Mabsout ME, Tassoulas LJ. Finite element model for the simulation of pile driving. Int J Numer Meth Eng 1994;37:257–78. [26] Helwany S. Applied Soil Mechanics with ABAQUS Applications. John Wiley & Sons; 2007. [27] Qiu G, Henke S, Grabe J. Application of a Coupled Eulerian-Lagrangian approach on geomechanical problems involving large deformations. Comput Geotech 2011;38(1):30–9. [28] Zhu SY, Fu Q, Cai CB. Damage evolution and dynamic response of cement asphalt mortar layer of slab track under vehicle dynamic load. Sci China Technol Sci 2014;57(10):1883–94. [29] Hu H-T, Huang C-S, Wu M-H, Wu Y-M. Nonlinear analysis of axially loaded concrete-filled tube columns with confinement effect. J Struct Eng 2003;129(10):1322–9. [30] Chai JC, Shen SL, Liu MD, Yuan DJ. Predicting the performance of embankments on PVD-improved subsoils. Comput Geotech 2018;93:222–31. [31] Yin ZY, Zhu QY, Zhu JG. Experimental investigation on creep behavior of soft clays: Review and development. Rock Soil Mech 2013;34(2):1–17. (in Chinese). [32] Augustesen A, Liingaard M, Lade Poul V. Evaluation of Time-Dependent Behavior of Soils. Int J Geomech 2004;4(3):137–56. [33] Almeida MSS, Maria PELS, Martins ISM, Spotti AP, Coelho LBM. Consolidation of a very soft clay with vertical drains. Géotechnique 2000;50(6):633–43. [34] Nogami T, li M. Consolidation of clay with a system of vertical and horizontal drains. J Geotech Geoenviron Eng 2003;129(9):838–48. [35] Indraratna B, Kan ME, Potts D, Rujikiatkamjorn C, Sloan SW. Analytical solution and numerical simulation of vacuum consolidation by vertical drains beneath circular embankments. Comput Geotech 2016;80:83–96. [36] Li Y, Chen YQ, Podlubny I. Stability of fractional-order nonlinear dynamic systems: Lyapunov direct method and generalized Mittag-Leffler stability. Comput Math Appl 2010;59(5):1810–21. [37] Zhang XW, Wang CM, Li JX. Experimental study of coupling behaviors of consolidation-creep of soft clay and its mechanism (in Chinese). Rock Soil Mech 2011;32(12):3584–90. [38] Yin ZY, Jin YF, Shen SL, Huang HW. An efficient optimization method for identifying parameters of soft structured clay by an enhanced genetic algorithm and elastic–viscoplastic model. 2017;12(4):849–67. [39] Furudoi T, Kobayashi M. Geotechnical issues and approach on Kansai International Airport project prediction and performance of settlement. Doboku Gakkai Ronbunshuu 2009;65(4):998-1017.

12