Constitutive Modeling of Physical Properties of Coastal Sand during Tunneling Construction Disturbance

: This paper presents a simple but workable constitutive model for the stress–strain relationship of sandy soil during the process of tunneling construction disturbance in coastal cities. The model was developed by linking the parameter K and internal angle ϕ of the Duncan–Chang model with the disturbed degree of sand, in which the effects of the initial void ratio on the strength deformation property of sands are considered using a uniﬁed disturbance function based on disturbed state concept theory. Three cases were analyzed to investigate the validity of the proposed constitutive model considering disturbance. After validation, the proposed constitutive model was further incorporated into a 3D ﬁnite element framework to predict the soil deformation caused by shield construction. It was found that the simulated results agreed well with the analytical solution, indicating that the developed numerical model with proposed constitutive relationship is capable of characterizing the mechanical properties of sand under tunneling construction disturbance.


Introduction
There is an ongoing demand for tunnel construction in coastal cities that need to shore up their crumbling infrastructure, seeking more efficient and less polluting modes of transportation. However, significant release of stress involved in the tunnel construction may cause catastrophic consequences for neighboring structures and underground works due to the excessive settlement and instability of the load-bearing soil layers. As reported by Chen et al. [1], the common lining uplift of Ningbo Metro Line 1 in eastern China during the tunneling construction stage reached more than 30 mm, which resulted in local cracks in tunnel linings and surrounding buildings. Consequently, estimation of potential ground movement during tunnel constructions is of great importance for civil engineers involved in the safe design of tunnels and their construction [2][3][4][5][6]. As a validated approach, finite element methods (FEMs) have been widely adopted to estimate the deformation characteristics of ground associated with complicated tunneling excavation [7][8][9][10]. To make predictions accurate, the essential features of soil behavior have to be reproduced by using suitable constitutive models with an FEM [11]. Addenbrooke et al. [12] developed a 2D FEM to investigate tunnel-induced ground movements in which the nonlinear behavior of soils is reproduced by adopting the Duncan-Chang model. Later, Zhang et al. [13] extended this framework to a 3D analysis of nailed soil structures under working loads. Mroueh and Shahrour [14] developed a full 3D finite element model to study the interaction between tunneling in soft soils and adjacent structures based on an elastic perfectly plastic constitutive relation with a Mohr-Coulomb criterion. Karakus and Fowell [15] utilized . Figure 1. Shield construction disturbance to the mechanical properties of coastal sand.

Laboratory Test
The samples used in the tests were composed of ISO standard sand provided by the Xiamen Company of China. The particle size distribution (PSD) for the tested samples was within the range between 0.25 mm and 1 mm, as shown in Figure 2. The sample physical In the present research, a series of triaxial compression tests were conducted for dry and saturated sands with different initial void ratios. The tested results were used to modify the disturbance function in terms of K and ϕ, utilizing DSC theory. Based on the proposed disturbance function, a modified Duncan-Chang model [26] taking into account construction disturbance was established. The developed model was applied to reproduce the physical properties of another kind of sand in a disturbed state to identify the model's validity and effectiveness. The validated framework was further incorporated into a 3D finite element model to predict the soil deformation caused by shield construction.

Laboratory Test
The samples used in the tests were composed of ISO standard sand provided by the Xiamen Company of China. The particle size distribution (PSD) for the tested samples was within the range between 0.25 mm and 1 mm, as shown in Figure 2. The sample physical parameters are listed in Table 1 [23], where G s is the specific gravity; e max and e min are the maximum and minimum void ratio, respectively; w is the water content; and C u and C c are the coefficients of uniformity and curvature, respectively.
. Figure 1. Shield construction disturbance to the mechanical properties of coastal sand.

Laboratory Test
The samples used in the tests were composed of ISO standard sand provided by the Xiamen Company of China. The particle size distribution (PSD) for the tested samples was within the range between 0.25 mm and 1 mm, as shown in Figure 2. The sample physical parameters are listed in Table 1 [23], where Gs is the specific gravity; emax and emin are the maximum and minimum void ratio, respectively; w is the water content; and Cu and Cc are the coefficients of uniformity and curvature, respectively.  All the tests were performed in the geotechnical laboratory of Zhejiang University in China. The triaxial device used in this experiment was the SJ-1A model, designed and manufactured by Guo Dian Nanjing Automation Company Limited, China. In this apparatus, a servohydraulic system is applied to control the cyclic vertical stress and frequency of loading, whereas an oil pressure type piston is applied to control the confining pressure, which varies from 0 to 2 MPa.  All the tests were performed in the geotechnical laboratory of Zhejiang University in China. The triaxial device used in this experiment was the SJ-1A model, designed and manufactured by Guo Dian Nanjing Automation Company Limited, China. In this apparatus, a servohydraulic system is applied to control the cyclic vertical stress and frequency of loading, whereas an oil pressure type piston is applied to control the confining pressure, which varies from 0 to 2 MPa.
The specimens were prepared with four different initial void ratios (e 0 ). The specific values for these specimens are listed in Table 2, in which m and m a are parameters defining the mass of the specimen and the mass of each layer, respectively. The specimens in the tests were prepared following the techniques of the SL237-1999 standard [27] and tested at three different confining levels (100, 200 and 300 kPa) and a fixed value of e 0 . A total of 24 standard undrained monotonic triaxial tests were conducted by Zhu et al. [23] at a strain rate of 0.808 mm/min for dry sand, with the failure criterion being controlled by the peak strength. The test results are shown in Figure 3, in which σ 1 and σ 3 are the principal stresses corresponding to the axial and circumferential directions in the test, respectively, and ε 1z represents the axial strain. As can be seen, as the initial void ratio e 0 decreased, the stress-strain curve for any confining pressure became steeper and the peak strength increased. Moreover, the strain that corresponded to the peak strength for all tests was approximately within the range between 2% and 3%. The specimens were prepared with four different initial void ratios (e0). The specific values for these specimens are listed in Table 2, in which m and ma are parameters defining the mass of the specimen and the mass of each layer, respectively. The specimens in the tests were prepared following the techniques of the SL237-1999 standard [27] and tested at three different confining levels (100, 200 and 300 kPa) and a fixed value of e0. A total of 24 standard undrained monotonic triaxial tests were conducted by Zhu et al. [23] at a strain rate of 0.808 mm/min for dry sand, with the failure criterion being controlled by the peak strength. The test results are shown in Figure 3, in which σ1 and σ3 are the principal stresses corresponding to the axial and circumferential directions in the test, respectively, and ε1z represents the axial strain. As can be seen, as the initial void ratio e0 decreased, the stress-strain curve for any confining pressure became steeper and the peak strength increased. Moreover, the strain that corresponded to the peak strength for all tests was approximately within the range between 2% and 3%.

Initial Tangent Modulus and Internal Friction Angle for Different Void Ratios
Kondner et al. [28,29] suggested that the nonlinear behavior of soils such as clay and sand can be effectively estimated with a hyperbola function, expressed as where a and b are model parameters for determining the initial tangent modulus and critical stress state when the stress-strain curve approaches infinite strain (σ1 − σ3)ult. Duncan and Chang [26] suggested that a and b could be determined as where the subscripts 95% and 70% represent the ratio between the values of stress difference (σ1 − σ3) and their peak values of strength (σ1 − σ3)f, respectively.
The relationship between (σ1 − σ3)f and (σ1 − σ3)ult is established by the failure ratio Rf as

Initial Tangent Modulus and Internal Friction Angle for Different Void Ratios
Kondner et al. [28,29] suggested that the nonlinear behavior of soils such as clay and sand can be effectively estimated with a hyperbola function, expressed as where a and b are model parameters for determining the initial tangent modulus and critical stress state when the stress-strain curve approaches infinite strain (σ 1 − σ 3 ) ult . Duncan and Chang [26] suggested that a and b could be determined as where the subscripts 95% and 70% represent the ratio between the values of stress difference (σ 1 − σ 3 ) and their peak values of strength (σ 1 − σ 3 ) f , respectively. The relationship between (σ 1 − σ 3 ) f and (σ 1 − σ 3 ) ult is established by the failure ratio R f as With the test results in Figure 3, the parameters (σ 1 − σ 3 ) f , a, b, E i , (σ 1 − σ 3 ) ult and instant failure ratio R fi can be determined accordingly by Equations (2)-(4). The results are shown in Table 3. Based on the laboratory tests, Janbu [30] suggested that the relationship between the initial tangent modulus and the confining pressure could be expressed as where p a is the atmospheric pressure, K is a model constant and n is a dimensionless parameter related to the rate of variation of E i and σ 3 . Another form of Equation (5) can be stated as With regard to the Mohr-Coulomb failure criterion, the peak failure strength can be derived as [26][27][28][29][30][31] where c and ϕ are the cohesion and friction angle of the soil, respectively. Generally, c = 0 for sand, so the value ϕ can be derived as Then, the values of the parameters K, n and the instant internal friction angle ϕ i can be determined based on the experimental results in combination with Equations (6) and (8). The results are shown in Table 4, where ϕ and R f are the mean values of ϕ i and R fi associated with each confining pressure.
Basically, the parameters n and R f are not sensitive to the variation of e 0 defining the tunnel disturbance during the test. Therefore, only the relationships between K, ϕ and e 0 are provided, as shown in Figures 4 and 5, respectively. The test results show that both lnK and sinϕ are linearly proportional to the change of e 0 . Hence the relationship between K, ϕ and e 0 can be specified as follows: where d, f, g and h are dimensionless parameters which can be effectively determined using linear regression of the test results with correlation coefficients greater than 0.92. Basically, the parameters n and Rf are not sensitive to the variation of e0 defining the tunnel disturbance during the test. Therefore, only the relationships between K, φ and e0 are provided, as shown in Figures 4 and 5, respectively. The test results show that both lnK and sinφ are linearly proportional to the change of e0. Hence the relationship between K, φ and e0 can be specified as follows: where d, f, g and h are dimensionless parameters which can be effectively determined using linear regression of the test results with correlation coefficients greater than 0.92.
The parameter (σ1 − σ3)f defining the critical stress state of soil should also be modified by substituting Equation (10) into Equation (7) Basically, the parameters n and Rf are not sensitive to the variation of e0 defining the tunnel disturbance during the test. Therefore, only the relationships between K, φ and e0 are provided, as shown in Figures 4 and 5, respectively. The test results show that both lnK and sinφ are linearly proportional to the change of e0. Hence the relationship between K, φ and e0 can be specified as follows: where d, f, g and h are dimensionless parameters which can be effectively determined using linear regression of the test results with correlation coefficients greater than 0.92.
The parameter (σ1 − σ3)f defining the critical stress state of soil should also be modified by substituting Equation (10) into Equation (7)  E i can be determined by substituting Equation (9) into Equation (5) as The parameter (σ 1 − σ 3 ) f defining the critical stress state of soil should also be modified by substituting Equation (10) into Equation (7), with a simplified expression, as

Unified Disturbance Function
In this study, a generalized disturbance function D varying from −1 to 1 was used to determine the degree of soil disturbance during the tunnel construction, as shown in Figure 6. Basically, there is no disturbance in soil when it is at the initial void ratio state. As the soil becomes looser, the void ratio e increases, whereas the corresponding degree of disturbance (D) decreases. When e approaches e max , D approaches −1 and the sand arrives at the loosest state. In turn, as the sand becomes denser, the amount of e decreases, whereas D increases. When e approaches e min , D reaches 1 and the backfill arrives at the densest state. The relationship between the disturbed degree and void ratio can be expressed as below.

Unified Disturbance Function
In this study, a generalized disturbance function D varying from −1 to 1 was used to determine the degree of soil disturbance during the tunnel construction, as shown in Figure 6. Basically, there is no disturbance in soil when it is at the initial void ratio state. As the soil becomes looser, the void ratio e increases, whereas the corresponding degree of disturbance (D) decreases. When e approaches emax, D approaches −1 and the sand arrives at the loosest state. In turn, as the sand becomes denser, the amount of e decreases, whereas D increases. When e approaches emin, D reaches 1 and the backfill arrives at the densest state. The relationship between the disturbed degree and void ratio can be expressed as below.

Simplified Constitutive Model Considering Disturbance
(1) For "positive disturbance" (e ≤ e 0 ), Equation (13a) can be rewritten as Substituting Equation (14) into Equation (9), the parameter K D during the disturbance can be expressed as Substituting Equation (14) into Equation (10), the parameter ϕ D during disturbance can be expressed as Substituting Equations (4), (15) and (16) into Equation (1), the relationship between the deviator stress and axial strain considering "positive disturbance" can be expressed as (2) For "negative disturbance" (e > e 0 ), Equation (13b) can be rewritten as Substituting Equation (18) into Equation (9), the parameter K D during the disturbance can be expressed as Substituting Equation (18) into Equation (10), the parameter ϕ D during the disturbance can be expressed as Substituting Equations (4), (19) and (20) into Equation (1), the relationship between deviator stress and axial strain considering "negative disturbance" can be expressed as A total of twelve parameters, namely e 0 , e max , e min , R f0 , K 0 , n 0 , ϕ 0 , σ 3 , d, f, g and h, are covered in the proposed model. Among these, e 0 , e max and e min can be easily determined by the fundamental physical test of sand; ϕ 0 , k 0 , R f0 , n 0 , are the same as in the original Duncan-Chang model; and d, f, g and h can be calibrated by the traditional undrained triaxial tests of sand.

Verification
We next took another kind of dry sand, Fujian standard sand (FJ sand), as the test material and conducted a series of triaxial compression tests to assess the validity of the proposed simplified constitutive model for disturbed states. The particle size distributions of the Fujian standard sand mainly ranged from 0.25 mm to 1 mm, as shown in Figure 7 [23]. Its physical parameters are listed in Table 5. The associated parameters of the proposed simplified constitutive model are listed in Table 6. Suppose that the initial void ratio e 0 of the sand is 0.76 and the sands at other void ratios (e.g., e = 0.79, 0.73, 0.70) are at different disturbed states, with their corresponding disturbed degrees as shown in Table 7. The predicted results of the proposed model are shown in Figure 8a-d. It can be seen that the stress-strain relationship of the sandy soil was significantly affected, either positively or negatively, by the disturbances. The predicted results always agreed well with the test curve at any disturbed state.
Substituting Equation (18) into Equation (9), the parameter KD during the disturbance can be expressed as Substituting Equation (18) into Equation (10), the parameter φD during the disturbance can be expressed as Substituting Equations (4), (19) and (20) into Equation (1), the relationship between deviator stress and axial strain considering "negative disturbance" can be expressed as A total of twelve parameters, namely e0, emax, emin, f0 R , K0, n0, φ0, σ3, d, f, g and h, are covered in the proposed model. Among these, e0, emax and emin can be easily determined by the fundamental physical test of sand; φ0, k0, f0 R , n0, are the same as in the original Duncan-Chang model; and d, f, g and h can be calibrated by the traditional undrained triaxial tests of sand.

Verification
We next took another kind of dry sand, Fujian standard sand (FJ sand), as the test material and conducted a series of triaxial compression tests to assess the validity of the proposed simplified constitutive model for disturbed states. The particle size distributions of the Fujian standard sand mainly ranged from 0.25 mm to 1 mm, as shown in Figure 7 [23]. Its physical parameters are listed in Table 5. The associated parameters of the proposed simplified constitutive model are listed in Table 6. Suppose that the initial void ratio e0 of the sand is 0.76 and the sands at other void ratios (e.g., e = 0.79, 0.73, 0.70) are at different disturbed states, with their corresponding disturbed degrees as shown in Table  7. The predicted results of the proposed model are shown in Figure 8a-d. It can be seen that the stress-strain relationship of the sandy soil was significantly affected, either positively or negatively, by the disturbances. The predicted results always agreed well with the test curve at any disturbed state.

Application
The developed constitutive framework was further incorporated into the commercially available software ABAQUS to reproduce the ground movement during tunnel constructions and compare the acquired simulation with the analytical solution.
In this study, a metro tunnel was considered as being 26.2 km long, 6.2 m in diameter and 0.35 m thick. The computing parameters of the tunnel are listed in Table 8. Assuming that the soil was of isotropic and homogeneous behavior, only half of the whole tunnel (Figure 9) was modeled to optimize the computational cost. The FEM mesh consisted of 10,010 nodes and 9480 elements (six-and eight-node linear brick), conditioned with appropriate boundary conditions. The upper surface corresponding to the effective ground was free to move, for which the pressure induced by the self-gravity of the EPB-S machine

Application
The developed constitutive framework was further incorporated into the commercially available software ABAQUS to reproduce the ground movement during tunnel constructions and compare the acquired simulation with the analytical solution.
In this study, a metro tunnel was considered as being 26.2 km long, 6.2 m in diameter and 0.35 m thick. The computing parameters of the tunnel are listed in Table 8. Assuming that the soil was of isotropic and homogeneous behavior, only half of the whole tunnel ( Figure 9) was modeled to optimize the computational cost. The FEM mesh consisted of 10,010 nodes and 9480 elements (six-and eight-node linear brick), conditioned with appropriate boundary conditions. The upper surface corresponding to the effective ground was free to move, for which the pressure induced by the self-gravity of the EPB-S machine was taken into account by applying a constant distributed load equal to 20 kPa [1,32] (additional thrust p). The interactive behavior of the shield-soil wall due to the effects of fluid injections from the shield head was conditioned with an additional friction force τ (45 kPa) [33,34]. The physical and mechanical parameters of the soil for the analytical solution [23] are reported in Table 9.  [33,34]. The physical and mechanical parameters of the soil for the analytical solution [23] are reported in Table 9.   The soil was considered to be isotopically nonlinearly elastic and homogeneous, with Poisson's ratio constant during the shield construction disturbance. The additional thrust p and additional friction force τ were also assumed to be acting straight on the soil in contrast with the analytical solution [23]. The physical properties assumed for the soils are reported in Table 10.  Figures 10 and 11 show the stress and vertical displacement fields before additional thrust and additional friction force act upon the soil. The initial stress field is well balanced because the stress of the soil is in line with the depth and the maximum vertical displacement is 10 −7 mm.  The soil was considered to be isotopically nonlinearly elastic and homogeneous, with Poisson's ratio constant during the shield construction disturbance. The additional thrust p and additional friction force τ were also assumed to be acting straight on the soil in contrast with the analytical solution [23]. The physical properties assumed for the soils are reported in Table 10.  10 and 11 show the stress and vertical displacement fields before additional thrust and additional friction force act upon the soil. The initial stress field is well balanced because the stress of the soil is in line with the depth and the maximum vertical displacement is 10 −7 mm.   Figures 12 and 13 show the stress and vertical displacement fields after the additional thrust acts on the system. The vertical displacements before and behind the shield area show eminence and subsidence, respectively. However, there is no significant change in the stress field of the soil.    Figures 12 and 13 show the stress and vertical displacement fields after the additional thrust acts on the system. The vertical displacements before and behind the shield area show eminence and subsidence, respectively. However, there is no significant change in the stress field of the soil.   Figures 12 and 13 show the stress and vertical displacement fields after the additional thrust acts on the system. The vertical displacements before and behind the shield area show eminence and subsidence, respectively. However, there is no significant change in the stress field of the soil.   Figures 12 and 13 show the stress and vertical displacement fields after the additional thrust acts on the system. The vertical displacements before and behind the shield area show eminence and subsidence, respectively. However, there is no significant change in the stress field of the soil.    Figures 14 and 15 show the stress and vertical displacement fields after the additional friction force acts on the system. The vertical displacements before and behind the shield area also show eminence and subsidence, respectively. Moreover, there is significant change in the stress field of the soil.  During shield construction, a uniformed change in the relative density of the soil is often considered. Different states of soil around the tunnel during shield construction can then be covered by assuming five different relative soil densities, in which Dr = 0.5 is the initial state, Dr = 0.3 and 0.4 are the negative disturbed states, and Dr = 0.6 and 0.7 are the  Figures 14 and 15 show the stress and vertical displacement fields after the additional friction force acts on the system. The vertical displacements before and behind the shield area also show eminence and subsidence, respectively. Moreover, there is significant change in the stress field of the soil.      During shield construction, a uniformed change in the relative density of the soil is often considered. Different states of soil around the tunnel during shield construction can then be covered by assuming five different relative soil densities, in which D r = 0.5 is the initial state, D r = 0.3 and 0.4 are the negative disturbed states, and D r = 0.6 and 0.7 are the positive disturbed states. The corresponding degree of disturbance computed by Equation (14) for each case is shown in Table 11. The additional thrust and friction force at the different disturbed states can be calculated using the 3D finite element model, taking into account ground movements and considering disturbance. As a comparison, the results based on the analytical solution proposed by Zhu et al. [23] are incorporated in Figures 16 and 17. As can be seen, the predicted vertical ground movement of the 3D FEM always agreed well with the analytical solution in any disturbed state. The observations in Figures 10-17 indicate that the proposed constitutive model can characterize the mechanical properties of sand under construction disturbance. positive disturbed states. The corresponding degree of disturbance computed by Equation (14) for each case is shown in Table 11. The additional thrust and friction force at the different disturbed states can be calculated using the 3D finite element model, taking into account ground movements and considering disturbance. As a comparison, the results based on the analytical solution proposed by Zhu et al. [23] are incorporated in Figures 16  and 17. As can be seen, the predicted vertical ground movement of the 3D FEM always agreed well with the analytical solution in any disturbed state. The observations in Figures  10-17 indicate that the proposed constitutive model can characterize the mechanical properties of sand under construction disturbance.

Conclusions
The motivation of the present study was to develop a practical constitutive model with the capacity to predict nonlinear soil behavior during the tunneling construction disturbance process. First, a series of undrained triaxial tests were conducted for samples of ISO standard sand with different initial void ratios. It was found that both the slope of the stress-strain curve and the peak strength increase when the value of e0 decreases. Based on the test results, a unified disturbance function was proposed based on DSC theory, in which the void ratio was selected as the disturbance parameter. Then, a novel approach relating the parameter K and internal angle φ to the disturbed degree was derived to modify the constitutive model considering construction disturbance. The proposed constitutive model was used to predict the physical properties of Fujian standard sand in a disturbed state. The results show that the predicted results for the stress-strain relationship of the soil always agreed well with the experimental data for any disturbed state. Finally, the proposed constitutive model was incorporated into the finite element modeling and validated against the analytical solution [23]. The results show that the predicted results in terms of vertical ground movements were in good agreement with the analytical solution [23] for any disturbed state, indicating that the developed model is capable of reproducing the mechanical behavior of sandy soil across the whole process of shield construction and can be extensively implemented for predicting construction disturbance effects in practical tunneling engineering.
It should be noted here that the proposed disturbance function is fundamental and does not cover physical parameters such as water content, mass density, etc. Moreover,

Conclusions
The motivation of the present study was to develop a practical constitutive model with the capacity to predict nonlinear soil behavior during the tunneling construction disturbance process. First, a series of undrained triaxial tests were conducted for samples of ISO standard sand with different initial void ratios. It was found that both the slope of the stress-strain curve and the peak strength increase when the value of e 0 decreases. Based on the test results, a unified disturbance function was proposed based on DSC theory, in which the void ratio was selected as the disturbance parameter. Then, a novel approach relating the parameter K and internal angle ϕ to the disturbed degree was derived to modify the constitutive model considering construction disturbance. The proposed constitutive model was used to predict the physical properties of Fujian standard sand in a disturbed state. The results show that the predicted results for the stress-strain relationship of the soil always agreed well with the experimental data for any disturbed state. Finally, the proposed constitutive model was incorporated into the finite element modeling and validated against the analytical solution [23]. The results show that the predicted results in terms of vertical ground movements were in good agreement with the analytical solution [23] for any disturbed state, indicating that the developed model is capable of reproducing the mechanical behavior of sandy soil across the whole process of shield construction and can be extensively implemented for predicting construction disturbance effects in practical tunneling engineering.
It should be noted here that the proposed disturbance function is fundamental and does not cover physical parameters such as water content, mass density, etc. Moreover, the inherent relationship between the mechanical parameters, such as stiffness, cohesion and internal friction angle, of clay or sand-clay admixture and the aforementioned physical parameters under the disturbed state of tunneling construction should be further addressed. More advanced solutions including these parameters will be further studied in future research.