Influence of Shield Attitude Change on Shield – Soil Interaction

The mechanism of shield–soil interaction and multi-phase equilibrium control theory in shield tunneling process still lack sufficient understanding. Aiming at this problem, with the improved calculation method of loose earth pressure, the initial boundary problem of shield attitude calculation was solved. Based on the ground reaction curve, the shield–soil interaction was simulated by the equivalent springs, and the displacement of surrounding soil was calculated during the change of the shield attitude. Then, the theoretical method of surrounding soil load acting on the shield were obtained. In summary, the calculation method of shield attitude was obtained. This method has three main applications in engineering, namely the inversion of shield–soil interaction force, the prediction of pitch angle and the prediction of yawing angle. Finally, combined with Jinan Metro Line R2 shield tunnel project, the shield attitude was monitored in real time and compared with the theoretical value. The results show that the trend of the theoretical values of pitch angle and yawing angle were basically the same as the measured value, but the theoretical value was generally larger than the measured value. The research results can provide a useful reference for the shield attitude adjustment.


Introduction
Shield tunneling has been widely used in the construction of urban subway tunnels.Due to the uncertainty of load and stratum during the excavation process, the shield cannot always be consistent with the theoretical design axis [1][2][3][4].The automatic control system used to control the shield attitude is based on existing engineering experience and not supported by appropriate theories [5].Therefore, how to properly control the advance load is one of the main concerns in the construction of shield tunnels.At present, relevant research mainly focuses on the two directions of advance load modeling and automatic correction.
The difficulty in establishing a shield tunneling mechanics model is that the forces around the shield cannot be calculated and expressed accurately.Therefore, many scholars have carried out a series of studies on the establishment of shield tunneling mechanics models and attitude control.Koyama et al. [6] mentioned that it is impossible to deal with the precise control of the shield machine in the case of complex geological structures and sharp curves.Xie et al. [7] applied the automatic trajectory tracking control system to the control of the shield's route, but these systems also only rely on empirical relationships, without precise theoretical background.Shimizu et al. [8] analyzed the tunneling and motion laws of tunnel construction without the influence of nonlinear factors.Shimizu et al. [9] established a mathematical analysis model of shield tunneling movement.Robust theory was applied to shield attitude control by Komiya et al. [10].Based on the force balance condition of the shield, Sugimoto et al. [11] determined the relationship between the advance load and the earth pressure, and then connected the earth pressure with the shield attitude, thus establishing the theoretical model of the shield advance load.Yue et al. [12] and Sun et al. [13] proposed a shield attitude and trajectory automatic control system based on the sliding mode robust controller and a shield attitude dynamic coordination control system based on the dynamic model of the shield attitude adjustment process.Ates et al. [14] believed that it is crucially important to select a proper TBM and define its basic specifications such as installed cutter-head torque and TBM thrust capacities.Zhang et al. [15] established a predicting model of the thrust and torque for the total load that fully reflects the influences of geological, operating, and structural parameters.In summary, the current research on shield attitude is more about the application of control theory, while ignoring the basic problems such as shield-soil interaction and shield tunneling mechanics.
With the necessity to accurately predict performance of machines in different ground conditions, many researchers have worked to develop new prediction models or adjustment factors for the common existing models [16].Acaroglu et al. [17] established a model to predict specific energy requirement of constant cross-section disc cutters in the rock cutting process by using fuzzy logic method.Barton et al. [18] proposed a model named Q TBM that uses many input parameters (such as RQD, joint condition, Stress condition, intact rock strength, quartz content and TBM thrust).The modified CSM model was added rock mass properties as input parameters into the model [19].In hard rock formations, such models have often shown results that it would not expand too much.However, in the soft soil layer, the interaction model for shield tunneling is still in its infancy.
In view of the above deficiencies, the influence of the shield attitude parameter on the shield-soil interaction was studied.The shield-soil interaction was represented by an equivalent spring, and the interaction force between the shield and the soil was solved and analyzed according to the ground reaction curve.Based on this, the theoretical calculation method of shield attitude was proposed and compared with the actual project.

Model of Load Acting on Shield
Based on the shield load model established by Sugimoto [20], the contact relationship between the shield and the stratum is blurred.The force of shield-soil interaction is decomposed along the axis direction, the normal direction, and the tangential direction of the shield, and the load model is shown in Figure 1.
The loads acting on the shield are: f 1 , shield self-weight; f 2 , force on the shield tail; f 3 , advance load provided by the jacks; f 4 , load acting on the cutter head; and f 5 , surrounding soil load acting on the shield periphery [21].
Appl.Sci.2019, 9, x FOR PEER REVIEW 2 of 20 et al. [9] established a mathematical analysis model of shield tunneling movement.Robust theory was applied to shield attitude control by Komiya et al. [10].Based on the force balance condition of the shield, Sugimoto et al. [11] determined the relationship between the advance load and the earth pressure, and then connected the earth pressure with the shield attitude, thus establishing the theoretical model of the shield advance load.Yue et al. [12] and Sun et al. [13] proposed a shield attitude and trajectory automatic control system based on the sliding mode robust controller and a shield attitude dynamic coordination control system based on the dynamic model of the shield attitude adjustment process.Ates et al. [14] believed that it is crucially important to select a proper TBM and define its basic specifications such as installed cutter-head torque and TBM thrust capacities.Zhang et al. [15] established a predicting model of the thrust and torque for the total load that fully reflects the influences of geological, operating, and structural parameters.In summary, the current research on shield attitude is more about the application of control theory, while ignoring the basic problems such as shield-soil interaction and shield tunneling mechanics.
With the necessity to accurately predict performance of machines in different ground conditions, many researchers have worked to develop new prediction models or adjustment factors for the common existing models [16].Acaroglu et al. [17] established a model to predict specific energy requirement of constant cross-section disc cutters in the rock cutting process by using fuzzy logic method.Barton et al. [18] proposed a model named QTBM that uses many input parameters (such as RQD, joint condition, Stress condition, intact rock strength, quartz content and TBM thrust).The modified CSM model was added rock mass properties as input parameters into the model [19].In hard rock formations, such models have often shown results that it would not expand too much.However, in the soft soil layer, the interaction model for shield tunneling is still in its infancy.
In view of the above deficiencies, the influence of the shield attitude parameter on the shieldsoil interaction was studied.The shield-soil interaction was represented by an equivalent spring, and the interaction force between the shield and the soil was solved and analyzed according to the ground reaction curve.Based on this, the theoretical calculation method of shield attitude was proposed and compared with the actual project.

Model of Load Acting on Shield
Based on the shield load model established by Sugimoto [20], the contact relationship between the shield and the stratum is blurred.The force of shield-soil interaction is decomposed along the axis direction, the normal direction, and the tangential direction of the shield, and the load model is shown in Figure 1.
The loads acting on the shield are: f1, shield self-weight; f2, force on the shield tail; f3, advance load provided by the jacks; f4, load acting on the cutter head; and f5, surrounding soil load acting on the shield periphery [21].If the shield is regarded as a rigid body, the distance between any two points on the shield remains unchanged during the tunneling process.Six parameters are required to describe the shield attitude: cutterhead center coordinate (x 0 , y 0 , z 0 ), yawing angle α, pitch angle β, and rolling angle Ω.The attitude angles are shown in Figure 2a.
The coordinate system is established as shown in Figure 2b: in the geodetic coordinate system C E , the x-axis is vertically downward and the y-axis is on the horizontal plane; in the shield-structure coordinate system C M , the p-axis is vertically downward and the r-axis is along the shield axis; and the coordinate system C M is rotated by a certain angle around the r-axis to obtain a coordinate system C MR .
Appl.Sci.2019, 9, x FOR PEER REVIEW 3 of 20 If the shield is regarded as a rigid body, the distance between any two points on the shield remains unchanged during the tunneling process.Six parameters are required to describe the shield attitude: cutterhead center coordinate (x0, y0, z0), yawing angle α, pitch angle β, and rolling angle Ω.The attitude angles are shown in Figure 2a.
The coordinate system is established as shown in Figure 2b: in the geodetic coordinate system C E , the x-axis is vertically downward and the y-axis is on the horizontal plane; in the shield-structure coordinate system C M , the p-axis is vertically downward and the r-axis is along the shield axis; and the coordinate system C M is rotated by a certain angle around the r-axis to obtain a coordinate system C MR .
(a) (b) When the shield is considered to be in a static equilibrium state during the tunneling process, the static balance equation must be met: where e = p, q, and r represent different directions of force; and b = E, M, and MR are expressed in different coordinate systems

Initial Earth Pressure Calculation
One of the keys to establishing a shield-soil interaction model during the change of shield attitude is the calculation of initial earth pressure.Loose earth pressure theory is often used for deformation analysis of segments.The initial shield stress state should also be in the non-uniform state.This part solves the initial boundary problem of the shield interaction model by means of relevant theories.

Calculation of loose soil pressure in the overlaying soil
Terzaghi deduced the classical theory of loose earth pressure through the trap-door test [22].The analytical expression of the vertical stress σv of the trap-door is: When the shield is considered to be in a static equilibrium state during the tunneling process, the static balance equation must be met: Then: where e = p, q, and r represent different directions of force; and b = E, M, and MR are expressed in different coordinate systems.

Initial Earth Pressure Calculation
One of the keys to establishing a shield-soil interaction model during the change of shield attitude is the calculation of initial earth pressure.Loose earth pressure theory is often used for deformation analysis of segments.The initial shield stress state should also be in the non-uniform state.This part solves the initial boundary problem of the shield interaction model by means of relevant theories.

Calculation of Loose Soil Pressure in the Overlaying Soil
Terzaghi deduced the classical theory of loose earth pressure through the trap-door test [22].The analytical expression of the vertical stress σ v of the trap-door is: where σ v is the loose earth pressure; K 0 is the later pressure coefficient; H is the overlying soil thickness; P 0 is the upper load; c is the cohesion; ϕ is the internal friction angle; γ is the weight density; and B 1 is the loose band width.
where R is the shield radius.
Li [23] introduced the factor of formation loss rate based on the stress of Terzaghi loose soil, and set the key parameter A 1 to correct the theoretical solution of Terzaghi loose soil pressure: (5) where K a is the active side pressure coefficient, K a = tan 2 (45 • -ϕ/2).
According to Figure 3, the following conclusions can be drawn: where s is the displacement of the dome.
Appl.Sci.2019, 9, x FOR PEER REVIEW 4 of 20 where σv is the loose earth pressure; K0 is the later pressure coefficient; H is the overlying soil thickness; P0 is the upper load; c is the cohesion; φ is the internal friction angle; γ is the weight density; and B1 is the loose band width.
where R is the shield radius.
Li [23] introduced the factor of formation loss rate based on the stress of Terzaghi loose soil, and set the key parameter A1 to correct the theoretical solution of Terzaghi loose soil pressure: ( ) where Ka is the active side pressure coefficient, Ka = tan 2 (45°-φ / 2).According to Figure 3, the following conclusions can be drawn: where s is the displacement of the dome.

Initial force acting on the shield periphery
The earth pressure around the shield increases with the increase of the depth of the soil, as shown in Figure 4. η is the angle with the p-axis in the C M coordinate system.To simplify the calculation, linear regression method is adopted.In the process of η from π/2 to 3π/2, the corresponding normal earth pressure acting on the shield at any point is linearly increased, so that the normal force on the shield can be obtained.The normal earth pressure at any point in the shell is shown in Figure 5.
Therefore, the earth pressure on the shield shell can be calculated as follows:

Initial Force Acting on the Shield Periphery
The earth pressure around the shield increases with the increase of the depth of the soil, as shown in Figure 4. η is the angle with the p-axis in the C M coordinate system.To simplify the calculation, linear regression method is adopted.In the process of η from π/2 to 3π/2, the corresponding normal earth pressure acting on the shield at any point is linearly increased, so that the normal force on the shield can be obtained.The normal earth pressure at any point in the shell is shown in Figure 5.
Therefore, the earth pressure on the shield shell can be calculated as follows: Appl.Sci.2019, 9, 1812 5 of 19 where p 1 is the normal earth pressure at the top of the shield, calculated according to improved loose soil pressure formula [23]; p 2 is the normal earth pressure at the bottom of the shield; D is the shield diameter; and L is the shield length.
where p1 is the normal earth pressure at the top of the shield, calculated according to improved loose soil pressure formula [23]; p2 is the normal earth pressure at the bottom of the shield; D is the shield diameter; and L is the shield length.

Basic assumption
The shield is constrained by the soil.Based on the characteristics of the ground reaction curve, the soil can be approximated as an equivalent spring within a certain range (Figure 6).The equivalent spring is nonlinear deformation, and the correlation coefficient can be obtained from the ground reaction curve [24].According to the deformation of the soil spring, the interaction force between the shield and the soil can be obtained.
The center of gravity of the shield machine is generally not in the geometric centroids, as shown in Figure 7a.Point S is the center of gravity of the shield, Point O is the geometric centroid, and ls is the distance between two points.
When the shield attitude is adjusted in the actual project, the yawing angle and pitch angle should be changed simultaneously.However, when the time interval of shield attitude adjustment tends to infinitesimal, it can be considered that the pitch angle and yawing angle are in succession.The order in which the pitch and yawing angles occur does not affect the final calculation result from where p1 is the normal earth pressure at the top of the shield, calculated according to improved loose soil pressure formula [23]; p2 is the normal earth pressure at the bottom of the shield; D is the shield diameter; and L is the shield length.

Basic assumption
The shield is constrained by the soil.Based on the characteristics of the ground reaction curve, the soil can be approximated as an equivalent spring within a certain range (Figure 6).The equivalent spring is nonlinear deformation, and the correlation coefficient can be obtained from the ground reaction curve [24].According to the deformation of the soil spring, the interaction force between the shield and the soil can be obtained.
The center of gravity of the shield machine is generally not in the geometric centroids, as shown in Figure 7a.Point S is the center of gravity of the shield, Point O is the geometric centroid, and ls is the distance between two points.
When the shield attitude is adjusted in the actual project, the yawing angle and pitch angle should be changed simultaneously.However, when the time interval of shield attitude adjustment tends to infinitesimal, it can be considered that the pitch angle and yawing angle are in succession.The order in which the pitch and yawing angles occur does not affect the final calculation result from

Basic Assumption
The shield is constrained by the soil.Based on the characteristics of the ground reaction curve, the soil can be approximated as an equivalent spring within a certain range (Figure 6).The equivalent spring is nonlinear deformation, and the correlation coefficient can be obtained from the ground reaction curve [24].According to the deformation of the soil spring, the interaction force between the shield and the soil can be obtained.
The center of gravity of the shield machine is generally not in the geometric centroids, as shown in Figure 7a.Point S is the center of gravity of the shield, Point O is the geometric centroid, and l s is the distance between two points.
When the shield attitude is adjusted in the actual project, the yawing angle and pitch angle should be changed simultaneously.However, when the time interval of shield attitude adjustment tends to infinitesimal, it can be considered that the pitch angle and yawing angle are in succession.The order in which the pitch and yawing angles occur does not affect the final calculation result from the aspect of the mechanics analysis.Assuming that the shield is performing attitude adjustment, it can be divided into four stages, as follows: (i) Stage 1: Gravity stage.The point of gravity action is moved to the geometric centroid while the gravity eccentric moment M G is generated, and only gravity is considered at this stage.Under the action of gravity, the shield machine produces a displacement of ∆s v in the vertical direction, and ∆s v is positive downward and negative upward, as shown in Figure 7b.
(ii) Stage 2: Gravity eccentric moment stage.The shield deflects due to gravity eccentric bending moment on the vertical plane.The pitch angle caused by the M G is called β G , as shown in Figure 7c.The angle is positive when its direction is counterclockwise, and the angle is negative when its direction is clockwise.(iii) Stage 3: Vertical bending moment stage.Upper and lower partition jacks produce deflection moments M β A on the vertical plane.The pitch angle caused by the M β A is called β A , as shown in Figure 7c.The angle is positive when its direction is counterclockwise, and the angle is negative when its direction is clockwise.(iv) Stage 4: Horizontal bending moment stage.Under the action of horizontal deflection bending moment M α A , α which is called yawing angle is generated on the horizontal plane, as shown in Figure 7d.α is positive when its direction is counterclockwise and negative when its direction is clockwise.
the aspect of the mechanics analysis.Assuming that the shield is performing attitude adjustment, it can be divided into four stages, as follows: (i) Stage 1: Gravity stage.The point of gravity action is moved to the geometric centroid while the gravity eccentric moment MG is generated, and only gravity is considered at this stage.
Under the action of gravity, the shield machine produces a displacement of Δsv in the vertical direction, and Δsv is positive downward and negative upward, as shown in Figure 7b.(ii) Stage 2: Gravity eccentric moment stage.The shield deflects due to gravity eccentric bending moment on the vertical plane.The pitch angle caused by the MG is called βG, as shown in Figure 7c.The angle is positive when its direction is counterclockwise, and the angle is negative when its direction is clockwise.(iii) Stage 3: Vertical bending moment stage.Upper and lower partition jacks produce deflection moments M A β on the vertical plane.The pitch angle caused by the M A β is called βA, as shown in Figure 7c.The angle is positive when its direction is counterclockwise, and the angle is negative when its direction is clockwise.(iv) Stage 4: Horizontal bending moment stage.Under the action of horizontal deflection bending moment M A α , α which is called yawing angle is generated on the horizontal plane, as shown in Figure 7d.α is positive when its direction is counterclockwise and negative when its direction is clockwise.

Geometric parameters
Assuming that the section of the shield at r = L has a displacement Δsv in the vertical direction and Δsh in the horizontal direction, the displacement occurring at each point on the shield is calculated.As shown in Figure 8, D is the shield diameter, r is the shield radius, L is the length of the shield, O is the initial position of the shield, O' is the final position of the shield, and η is the angle with the p-axis in the C M coordinate system.
As shown in Figure 8, the following inference can be made:  shown in Figure 7c.The angle is positive when its direction is counterclockwise, and the angle is negative when its direction is clockwise.(iv) Stage 4: Horizontal bending moment stage.Under the action of horizontal deflection bending moment M A α , α which is called yawing angle is generated on the horizontal plane, as shown in Figure 7d.α is positive when its direction is counterclockwise and negative when its direction is clockwise.

Geometric parameters
Assuming that the section of the shield at r = L has a displacement Δsv in the vertical direction and Δsh in the horizontal direction, the displacement occurring at each point on the shield is calculated.As shown in Figure 8, D is the shield diameter, r is the shield radius, L is the length of the shield, O is the initial position of the shield, O' is the final position of the shield, and η is the angle with the p-axis in the C M coordinate system.
As shown in Figure 8, the following inference can be made:

Geometric Parameters
Assuming that the section of the shield at r = L has a displacement ∆s v in the vertical direction and ∆s h in the horizontal direction, the displacement occurring at each point on the shield is calculated.As shown in Figure 8, D is the shield diameter, r is the shield radius, L is the length of the shield, O is the initial position of the shield, O' is the final position of the shield, and η is the angle with the p-axis in the C M coordinate system.
As shown in Figure 8, the following inference can be made: Polar coordinate conversion: x = ρ cos(−η); y = ρ sin(−η).Bring in Equation ( 4) to get: Appl.Sci.2019, 9, 1812 MM < 0 indicates that the shield is subjected to active earth pressure, and MM > 0 indicates that the shield is subjected to passive earth pressure.Therefore, the vertical and horizontal displacement of the shield: ' 0 MM < indicates that the shield is subjected to active earth pressure, and ' 0 MM > indicates that the shield is subjected to passive earth pressure.Therefore, the vertical and horizontal displacement of the shield:

Solution of vertical and horizontal forces of shield-soil interaction
Through the ground reaction curve (Figure 9), considering the change of the shield attitude, the vertical and horizontal interaction forces between the shield and the soil are deduced and analyzed under the premise of small deformation.
The ground reaction curve [24], which shows the relationship between the ground displacement and the earth pressure acting on the shield periphery as shown in Figure 9, can be represented by where K is the coefficient of earth pressure, which is defined as the earth pressure σ divided by the initial vertical earth pressure σv0 (Ki=σi/σv0); U is the ground displacement; a is the gradient of functions K(U), which represents the coefficient of subgrade reaction ai; subscripts v and h are the vertical horizontal directions, respectively; and subscripts o, min, and max are the initial, lower limit and upper limit of the coefficient of earth pressure, respectively.

Solution of Vertical and Horizontal Forces of Shield-Soil Interaction
Through the ground reaction curve (Figure 9), considering the change of the shield attitude, the vertical and horizontal interaction forces between the shield and the soil are deduced and analyzed under the premise of small deformation.The shield is discretized, divided into n parts along the circumferential direction, and divided into m parts along the longitudinal direction, as shown in Figure 10 Vertical and horizontal stress of each unit are obtained as follows: ( ) ( ) where σv0ij and σh0ij are the initial vertical and horizontal earth pressures, respectively, and the initial earth pressure is calculated as the loose earth pressure.Therefore, the conclusions can be obtained as follows: 2π L The ground reaction curve [24], which shows the relationship between the ground displacement and the earth pressure acting on the shield periphery as shown in Figure 9, can be represented by Appl.Sci.2019, 9, 1812 8 of 19 where K is the coefficient of earth pressure, which is defined as the earth pressure σ divided by the initial vertical earth pressure σ v0 (K i = σ i /σ v0 ); U is the ground displacement; a is the gradient of functions K(U), which represents the coefficient of subgrade reaction a i ; subscripts v and h are the vertical horizontal directions, respectively; and subscripts o, min, and max are the initial, lower limit and upper limit of the coefficient of earth pressure, respectively.The shield is discretized, divided into n parts along the circumferential direction, and divided into m parts along the longitudinal direction, as shown in Figure 10.
into m parts along the longitudinal direction, as shown in Figure 10 Vertical and horizontal stress of each unit are obtained as follows: ( ) ( ) where σv0ij and σh0ij are the initial vertical and horizontal earth pressures, respectively, and the initial earth pressure is calculated as the loose earth pressure.Therefore, the conclusions can be obtained as follows: where Aij is the area of each unit of the shield shell, and Aij= It can be known from Equations ( 1) and ( 2):

Solution of bending moment of shield-soil interaction
Assuming that the pitch angle of the shield is β (β=βG+βA) and the yawing angle is α, then the vertical and horizontal displacement of each element on the shield are determined as follows: Vertical and horizontal stress of each unit are obtained as follows: where σ v0ij and σ h0ij are the initial vertical and horizontal earth pressures, respectively, and the initial earth pressure is calculated as the loose earth pressure.Therefore, the conclusions can be obtained as follows: where A ij is the area of each unit of the shield shell, and It can be known from Equations ( 1) and ( 2):

Solution of Bending Moment of Shield-Soil Interaction
Assuming that the pitch angle of the shield is β (β = β G + β A ) and the yawing angle is α, then the vertical and horizontal displacement of each element on the shield are determined as follows: where ∆s v and ∆s h are the vertical and horizontal displacements of the elements on different sections, respectively.x M r is the value of the r axis in the C M coordinate system.

Solution Process of Yawing Angle and Pitch Angle
According to the above calculation theory, the calculation method of the pitch angle and yawing angle of the shield can be summarized.The specific calculation process is shown in Figure 11.
The first stage is to input the basic parameters such as soil parameters and shield parameters, and establish the shield-soil interaction model.The relationship curve between the vertical displacement and vertical force on the shield (RCDV) is obtained.Through the RCDV, the vertical displacement ∆s v under gravity is obtained.In the second stage, the ∆s v obtained in the first stage and the unknown pitch angle are set into Equation ( 18), and the relationship curve between the pitch angle and moment on the vertical plane (RCPV) is obtained, which is obtained under the action of gravity eccentric bending moment.Through the RCPV, the pitch angle β G generated by the M G and the β A generated by the M β A can be obtained.In the third stage, ∆s v , β G , and β A are brought into the calculation program, and the unknown yawing angle value is continuously given.The relationship curve between the yawing angle and horizontal bending moment applied by the shield (RCYH) is obtained.Through the RCYH, the M α A is calculated, then the yawing angle α is obtained.Finally, ∆s v , β G , β A , and α are output.
sections, respectively. is the value of the r axis in the C coordinate system.
It can be known from Equations ( 1) and ( 2):

Solution Process of Yawing Angle and Pitch Angle
According to the above calculation theory, the calculation method of the pitch angle and yawing angle of the shield can be summarized.The specific calculation process is shown in Figure 11.
The first stage is to input the basic parameters such as soil parameters and shield parameters, and establish the shield-soil interaction model.The relationship curve between the vertical displacement and vertical force on the shield (RCDV) is obtained.Through the RCDV, the vertical displacement △sv under gravity is obtained.In the second stage, the △sv obtained in the first stage and the unknown pitch angle are set into Equation ( 18), and the relationship curve between the pitch angle and moment on the vertical plane (RCPV) is obtained, which is obtained under the action of gravity eccentric bending moment.Through the RCPV, the pitch angle βG generated by the MG and the βA generated by the M A β can be obtained.In the third stage, △sv, βG, and βA are brought into the calculation program, and the unknown yawing angle value is continuously given.The relationship curve between the yawing angle and horizontal bending moment applied by the shield (RCYH) is obtained.Through the RCYH, the M A α is calculated, then the yawing angle α is obtained.Finally, △ sv, βG, βA, and α are output.

Calculation Example
Based on the above derivation results, the following set of typical parameters for calculation was selected.Suppose a shield tunnel is excavated in a uniform sand layer.The soil weight γ is 19 kN/m 3 , the cohesion c is 0 kPa, and the internal friction angle φ is 30°.The shield diameter D is 12 m, the

Calculation Example
Based on the above derivation results, the following set of typical parameters for calculation was selected.Suppose a shield tunnel is excavated in a uniform sand layer.The soil weight γ is 19 kN/m 3 , the cohesion c is 0 kPa, and the internal friction angle ϕ is 30 • .The shield diameter D is 12 m, the length L is 10.8 m, the shield mass G shield is 3000 t, and the thickness of covering soil is 24 m.The displacement of the dome is 5 mm, and the eccentricity of gravity l s is 3.5 m.M β A = 20 MN•m and M α A = 75 MN•m.The ground reaction curve coefficient is shown in Table 1 [11].
Through the Matlab programming calculation, the shield attitude parameters could be calculated.As shown in Figure 12a, the weight of the shield machine was 3000 t, and ∆s v was obtained (∆s v = 8.1 mm).In the second stage, the sum of M G and M β A is 125 MN•m, then the pitch angle β could be obtained (Figure 12b) (β = 1.16 • ).In the third stage, it was found that M α A = 75 MN•m, and then the yawing angle α could be obtained (Figure 12c) (α = 0.7 • ).
underneath important regions such as Beijing-Shanghai high-speed railway, Beijing-Fuzhou expressway and Yufu River.The advancement of shield machine usually encountered a silty clay soil and loess during the tunneling process in the studied section.The geological profile of the encountered stratum is displayed in Figure 14.The construction process of Ring 504 segments was selected for analysis and research.The soil characteristic parameters of Ring 504 are shown in Figure 15.The third and fourth layers are the soil parameters of the shield excavation section.expressway and Yufu River.The advancement of shield machine usually encountered a silty clay soil and loess during the tunneling process in the studied section.The geological profile of the encountered stratum is displayed in Figure 14.The construction process of Ring 504 segments was selected for analysis and research.The soil characteristic parameters of Ring 504 are shown in Figure 15.The third and fourth layers are the soil parameters of the shield excavation section.

Shield Details
An earth pressure balanced shield (EPB) machine of 6.68 m in diameter was used to excavate in Jinan Metro Line R2.The length of shield is 9 m, and its weight is about 500 t.Cutter opening ratio is 35%, as shown in Figure 16a.The shield machine specifications are summarized in Table 2.The shield

Shield Details
An earth pressure balanced shield (EPB) machine of 6.68 m in diameter was used to excavate in Jinan Metro Line R2.The length of shield is 9 m, and its weight is about 500 t.Cutter opening ratio is 35%, as shown in Figure 16a.The shield machine specifications are summarized in Table 2.The shield machine advance system uses 22 sets of cylinders, which are divided into four groups according to the position.The number of cylinders in Groups A and C is 4, and the number of cylinders in Groups B and D is 7.The cylinder section shown in Figure 16b.
35%, as shown in Figure 16a.The shield machine specifications are summarized in Table 2.The shield machine advance system uses 22 sets of cylinders, which are divided into four groups according to the position.The number of cylinders in Groups A and C is 4, and the number of cylinders in Groups B and D is 7.The cylinder section is shown in Figure 16b.

Engineering Application I: Inversion Calculation of Soil-Shield Interaction Force
The soil layer above the section of Ring 504 is composed of backfill and loess, which is a nonuniform clay stratum.According to the method of this paper and in [27], initial force acting on the shield periphery can be calculated in the non-uniform clay stratum.As shown in Figure 15, the section of Ring 504 is also a non-uniform layer.The loess layer and silty clay each account for about half of the entire excavation surface and can be approximated by 1:1.
Through the guiding system and data monitoring system of the EPB, when the shield had been advancing at Ring 504, partition cylinder thrust, pitch angle and yawing angle were recorded every 1 min.The time to excavate Ring 504 was about 50 min, but for various reasons, it was interrupted

Engineering Application I: Inversion Calculation of Soil-Shield Interaction Force
The soil layer above the section of Ring 504 is composed of backfill and loess, which is a non-uniform clay stratum.According to the method of this paper and in [27], initial force acting on the shield periphery can be calculated in the non-uniform clay stratum.As shown in Figure 15, the section of Ring 504 is also a non-uniform layer.The loess layer and silty clay each account for about half of the entire excavation surface and can be approximated by 1:1.
Through the guiding system and data monitoring system of the EPB, when the shield had been advancing at Ring 504, partition cylinder thrust, pitch angle and yawing angle were recorded every 1 min.The time to excavate Ring 504 was about 50 min, but for various reasons, it was interrupted during the tunneling process.Continuous monitoring data were more valuable and informative.The longest continuous boring time of the ring was about 15 min.To avoid the influence of shield start and stop, this study selected continuous monitoring data of 10 min in the middle of this time.The data for a period were selected for analysis and comparison.The measured data of the pitch angle and the yawing angle are shown in Figure 17.The yawing angle varied greatly, while the pitch angle remained basically stable.The pitch angle and the yawing angle had been known through the guiding system.According to the flow chart of the shield attitude calculation (Figure 11), the force of shield-soil interaction during the tunneling process can be inverted.The calculation results are shown in Figure 18. Figure 18 shows the distribution of the total normal earth pressure around the shield periphery.The height of graphs represents the length of shield, the width of graphs the shield periphery from 0 • to 360 • .Along the line at 180 • , it can be found that the normal pressure is gradually decreasing from the bottom to the top of graphs, indicating that the pitch angle of the shield is positive.The gradient of the normal pressure reduction reflects the magnitude of the pitch angle.The stress nephogram on the left or right of the line at 180 • shows the direction of the yawing angle of the shield.

Engineering Application II: Shield Pitch Angle Prediction
The second engineering application is the prediction of the pitch angle.Combined with the Jinan Metro Line R2, when the advance force applied by the four cylinder partitions is known, the attitude change of the shield can be predicted in advance.The soil mass Wsoil with the same volume as the shield is about 600 t, which is larger than the self-weight of the shield machine.Therefore, under the premise of not applying any external force, there will be a phenomenon of "upward floating", that is, Δsv < 0. The first step is to calculate Δsv (Δsv =-13 mm), according to the method in Figure 11.
The Figure 20 is a comparative analysis of the measured and predicted values of the shield pitch angle during the study period.It can be seen in the figure that the theoretical value is generally larger than the measured value, indicating that the moment actively applied by the shield through the jack cylinder cannot fully act on the shield itself.When the advance cylinder acts on the shield, the friction at the joint and the characteristics of the shield itself may cause the deflection moment to decrease, resulting in a smaller pitch angle actually being formed.The change of the pitch angle of the shield is small, and the fluctuation of the stress nephogram is mainly due to the change of the yawing angle.The stress nephogram can be used to preliminarily determine the degree and the direction of deflection on the horizontal plane of the shield.The greater is the yawing angle, the more severe is the distortion of the stress nephogram.For example, there was a sudden increase of the yawing angle between time = 4 min and time = 5 min, as can be seen from the stress nephogram.

Engineering Application II: Shield Pitch Angle Prediction
The second engineering application is the prediction of the pitch angle.Combined with the Jinan Metro Line R2, when the advance force applied by the four cylinder partitions is known, the attitude change of the shield can be predicted in advance.The soil mass W soil with the same volume as the shield is about 600 t, which is larger than the self-weight of the shield machine.Therefore, under the premise of not applying any external force, there will be a phenomenon of "upward floating", that is, ∆s v < 0. The first step is to calculate ∆s v (∆s v = −13 mm), according to the method in Figure 11. Figure 20 is a comparative analysis of the measured and predicted values of the shield pitch angle during the study period.It can be seen in the figure that the theoretical value is generally larger than the measured value, indicating that the moment actively applied by the shield through the jack cylinder cannot fully act on the shield itself.When the advance cylinder acts on the shield, the friction at the joint and the characteristics of the shield itself may cause the deflection moment to decrease, resulting in a smaller pitch angle actually being formed.

Engineering Application III: Shield Yawing Angle Prediction
The third engineering application is the prediction of the yawing angle.The Δsv is calculated in the first stage.The theoretical value of the pitch angle at each moment is shown in Figure 20.The known Δsv and β are brought into the calculation process of the third stage to obtain the RCYH (Figure 21).The M A actively applied by the EPB can be calculated by the thrust difference between Group B and Group D of the advance system.The change curve of the M A is shown in Figure 22.Then, the theoretical value of the yawing angle can be obtained.The comparison between the theoretical value and the measured value is shown in Figure 23.The trends of theoretical and measured values are both similar to the trend of the M A .The M A is obtained by measurement.It indicates that the M A applied by the shield can act more fully on the shield than M A β .
Comparing the measured values with the theoretical values in Figure 23, it can be found that the trend of the theoretical value is basically similar to the measured value, but the theoretical value is generally larger than the measured value.In the calculation of the theoretical value, the deflection moment value is the initial data, and the moment used for the shield correction is reduced due to mechanical loss, etc., so that the final yawing angle formed by the shield is small.

Engineering Application III: Shield Yawing Angle Prediction
The third engineering application is the prediction of the yawing angle.The Δsv is calculated in the first stage.The theoretical value of the pitch angle at each moment is shown in Figure 20.The known Δsv and β are brought into the calculation process of the third stage to obtain the RCYH (Figure 21).The M A actively applied by the EPB can be calculated by the thrust difference between Group B and Group D of the advance system.The change curve of the M A is shown in Figure 22.Then, the theoretical value of the yawing angle can be obtained.The comparison between the theoretical value and the measured value is shown in Figure 23.The trends of theoretical and measured values are both similar to the trend of the M A .The M A is obtained by measurement.It indicates that the M A applied by the shield can act more fully on the shield than M A β .
Comparing the measured values with the theoretical values in Figure 23, it can be found that the trend of the theoretical value is basically similar to the measured value, but the theoretical value is generally larger than the measured value.In the calculation of the theoretical value, the deflection moment value is the initial data, and the moment used for the shield correction is reduced due to mechanical loss, etc., so that the final yawing angle formed by the shield is small.

Engineering Application III: Shield Yawing Angle Prediction
The third engineering application is the prediction of the yawing angle.The ∆s v is calculated in the first stage.The theoretical value of the pitch angle at each moment is shown in Figure 20.The known ∆s v and β are brought into the calculation process of the third stage to obtain the RCYH (Figure 21).The M α A actively applied by the EPB can be calculated by the thrust difference between Group B and Group D of the advance system.The change curve of the M α A is shown in Figure 22.Then, the theoretical value of the yawing angle can be obtained.The comparison between the theoretical value and the measured value is shown in Figure 23 23, it can be found that the trend of the theoretical value is basically similar to the measured value, but the theoretical value is generally larger than the measured value.In the calculation of the theoretical value, the deflection moment value is the initial data, and the moment used for the shield correction is reduced due to mechanical loss, etc., so that the final yawing angle formed by the shield is small.

Conclusions
This study revealed the mechanism of the interaction between shield and soil during the change of shield attitude.Through theoretical analysis, programming calculation and field monitoring, the calculation method of the pitch angle and yawing angle of the shield was studied and verified.The main achievements of this study are outlined as follows: (1) Based on the ground reaction force curve, the interaction between the shield and the soil was simulated by the equivalent spring, and the theoretical calculation method of f 5 was obtained in the change of the shield attitude.(2) The improved calculation method of loose earth pressure solved the initial boundary problem of the shield attitude calculation.Combined with the theoretical calculation method of f 5 , the calculation process of the shield attitude was formed.(3) Based on the monitoring data of Jinan Metro Line R2, the interaction between the shield and the soil during the construction of the shield was inverted.Through the stress nephogram, the stress concentration area of the shield can be judged to guide the next step of attitude adjustment.(4) Based on the measured data of the project, the theoretical calculation method of the pitch angle and yawing angle of the shield was verified.The study found that the theoretical value was close to the measured value, but, since the moment generated by the advance cylinders cannot fully act on the shield itself, the theoretical values of the shield attitude angles were generally greater than the measured values.(5) The model of shield-soil interaction has some benefits for upcoming projects.Firstly, this model can guide the shield operator to perform shield attitude correction.Secondly, it is possible to initially obtain a position where the attitude control of any shield tunnel construction is difficult for upcoming projects.Acknowledgments: The research work described herein was funded by the National Basic Research Program of China.This financial support is gratefully acknowledged.

Conflicts of Interest:
The authors declare no conflict of interest.

Figure 1 .
Figure 1.Shield load model: (a) Section A-A; (b) Section B-B; and (c) Section C-C.Figure 1. Shield load model: (a) Section A-A; (b) Section B-B; and (c) Section C-C.

Figure 1 .
Figure 1.Shield load model: (a) Section A-A; (b) Section B-B; and (c) Section C-C.Figure 1. Shield load model: (a) Section A-A; (b) Section B-B; and (c) Section C-C.

Figure 13 .
Figure 13.Plan view of the district division and geology in Jinan [26].

Figure 13 .
Figure 13.Plan view of the district division and geology in Jinan [26].

Figure 13 .
Figure 13.Plan view of the district division and geology in Jinan [26].

Figure 16 .Table 2 .
Figure 16.Cutterhead and advance system of shield: (a) shield cutterhead; and (b) shield cylinder partition.Table 2. Summary of the main specifications for the EPB.Shield Type EPB External diameter (m) 6.68 The length of shield (m) 9.0 Shield weight (t) Approximately 500 Shield eccentricity (m) Approximately 0.5 Cutter opening ratio (%) 35 Maximum total thrust (kN) 40,860 Number of propulsion cylinders 22 Maximum advance speed (mm/min) 80
M A β actively applied by the EPB can be calculated by the thrust difference between Group A and Group C of the advance system.The MG is 2.5 MN•m, which can be calculated.The change curve of the M A β +MG can be obtained during construction, as shown in Figure 19.When M A β +MG and Δsv are known, the theoretical values of the shield pitch angle during the study period can be obtained.The theoretical values are calculated by the M A β +MG.Therefore, the trend of theoretical values should be similar to the trend of the M A β +MG.However, the measured values did not have this trend, indicating that the M A β +MG did not fully act on shield attitude adjustment.

The M β A
actively applied by the EPB can be calculated by the thrust difference between Group A and Group C of the advance system.The M G is 2.5 MN•m, which can be calculated.The change curve of the M β A + M G can be obtained during construction, as shown in Figure 19.When M β A + M G and ∆s v are known, the theoretical values of the shield pitch angle during the study period can be obtained.The theoretical values are calculated by the M β A +M G .Therefore, the trend of theoretical values should be similar to the trend of the M β A + M G .However, the measured values did not have this trend, indicating that the M β A + M G did not fully act on shield attitude adjustment.

Figure 19 .
Figure 19.The change graph of M A β +MG.

Figure 20 .
Figure 20.Contrast between the measured and the theoretical value of the shield pitch angle

Figure 19 . 20 Figure 19 .
Figure 19.The change graph of M β A + M G .

Figure 20 .
Figure 20.Contrast between the measured and the theoretical value of the shield pitch angle

Figure 20 .
Figure 20.Contrast between the measured and the theoretical value of the shield pitch angle.
. The trends of theoretical and measured values are both similar to the trend of the M α A .The M α A is obtained by measurement.It indicates that the M α A applied by the shield can act more fully on the shield than M β A .Comparing the measured values with the theoretical values in Figure

Figure 21 .
Figure 21.The relationship between the yawing angle and M A α .

Figure 22 .
Figure 22.The change graph of M A α .

Figure 23 .
Figure 23.Contrast between the measured and the theoretical value of the shield yawing angle.

Figure 21 . 20 Figure 21 .
Figure 21.The relationship between the yawing angle and M α A .

Figure 22 .
Figure 22.The change graph of M A α .

Figure 23 .
Figure 23.Contrast between the measured and the theoretical value of the shield yawing angle.

Figure 22 .
Figure 22.The change graph of M A α .

Figure 23 .
Figure 23.Contrast between the measured and the theoretical value of the shield yawing angle.

Figure 23 .
Figure 23.Contrast between the measured and the theoretical value of the shield yawing angle.

AuthorFunding:
Contributions: X.S., D.-J.Y. and D.-L.J. contributed equally to this work.The research work described herein was funded by the National Basic Research Program of China (973 Program: 2015CB057800) and Fundamental Research Funds for the Central Universities of China (NO.2017YJS151).

Table 2 .
Summary of the main specifications for the EPB.