2. Materials and Methods
2.1 Data Simulation of Drift Ice Impacting Diversion Tunnel
The control equations that describe the nonlinear motion, deformation and impact problem fall into two categories: Lagrange method and Euler method. This study adopts Lagrange method to solve contact-impact problem. In the Lagrange method, the quality of particle must remain unchanged during the particle moves with the object.
Particles move in the overall motion system must follow the basic flow equation which satisfies the energy, mass and momentum conservation governing equation. In this study, the basic equations of drift ice impacting tunnel lining are shown as follows.
2.1.1 Mass conservation equation
where [TeX:] $$\rho_{0}(X, 0)$$ is the medium density of the initial model [TeX:] $$(t=0) . J(X, t)$$ is the Jacobin determinate of the model.
2.1.2 Momentum conservation equation
where [TeX:] $$\sigma_{i j}$$ is the current architectonic stress tensor. [TeX:] $$b_{i}$$ is the acting force on the unit mass of drift ice model and tunnel model. [TeX:] $$\ddot{u}_{i}$$ is the particles acceleration of the drift ice and tunnel model which is calculated by Eq. (3):
where [TeX:] $$v_{i}(X, t)$$ denotes the instantaneous velocity of the particle in time t and the coordinate of the particle is X.
2.1.3 Energy conservation equation
where [TeX:] $$\dot{w}^{\text {int }}$$ is the derivation of internal energy of the drift ice and tunnel model with response to time. [TeX:] $$D_{i j}$$ is the strain rate tensor which is calculated by Eq. (5).
where [TeX:] $$v_{i}$$ is the instantaneous velocity (m/s) of the particle in direction x. [TeX:] $$v_{j}$$ is the instantaneous velocity (m/s) of the particle in direction y.
2.1.4 Boundary condition
The tunnel model current architectonic surface force boundary condition is
where, [TeX:] $$t_{i}$$ is the stress vector on surface element, [TeX:] $$n_{j}$$ is the current architectonic directed surface element, [TeX:] $$\sigma_{i j}$$ and is the stress vector of the current architectonic. The calculation requires the relative motion between drift ice and tunnel satisfies momentum conservation equation in the whole solution range. But in the practical engineering project, it is almost impossible to solute the results directly. In order to solve this problem, we introduced the numerical computation method which considers the weak form of the differential equation. It could get the momentum equation only by the inner product and further deduce the virtual power theory equations of the model. After the finite element discretization, the joint displacement equation of model is obtained.
Based on the weighted residual method in which virtual velocity is chosen as weighting coefficient, the weak form of the momentum equation is obtained. It is
where [TeX:] $$\delta v_{i}(X)$$ is the virtual velocity which is solved by Eq. (8).
where, [TeX:] $$R_{0}$$ and [TeX:] $$C^{0}$$ are the vector and damping matrix in time 0, respectively, [TeX:] $$A_{t}$$ is the boundary surface of the model. Based on subsection integral principle, the balance formula of the contact force can be expressed as follow:
In this way, all particles velocity, acceleration, virtual velocity and strain rate could be expressed as the follow equation set:
Transform all the equations in the equation set (10) to matrix form and substitute them into virtual power formula (9) to get Eq. (11).
where,[TeX:] $$\ddot{U}$$denotes the system acceleration array, [TeX:] $$f^{\text {int }}$$denotes the system internal force array,[TeX:] $$f^{e x t}$$ denotes the system external force array, and M denotes the system mass array irrelevant to time. Their expressions are (12) to (14).
Solve the above Eq. (14) to get the particle displacement [TeX:] $$\mu_{I}$$ in the current time, then further obtained structure strain and stress at that moment.
2.2 Finite Element Modelling
2.2.1 Case study
In the case study, we chose 20 m in the No. 3 Pandaoling tunnel (CH62+959.78--78+682.934) to analyze the impact of the drift ice on the tunnel by finite element modelling method. Pandaoling tunnel is in the Datong River to Qinwangchuan region water diversion project in western China is 15.723 km in length, longitudinal slope of 1/1000, design flow of Q = 29 [TeX:] $$m^{3} / s$$ and increasing flow of Q = 34 [TeX:] $$m^{3} / s$$. It is pressureless diversion tunnel with mixed type lining. The first lining is made up of anchor arm, spay concrete, reinforced mesh and steel arch. The secondary lining is cash-in-place concrete and reinforced concrete. The tunnel with 4.2 m clear width, 4.4 m clear height has arch walls and sections with antiarch baseboard. It has semicircle roof with 2.1 m radius and antiarch baseboard with 9.75 m radium. The fillet is added at the sidewall and antiarch baseboard connection, it is 0.404 m in height and 0.337m in length. The early tunnel support tunnel section is shown in Fig. 1. The tunnel model is shown in Fig. 2.
Tunnel section primary support.
2.2.2 Project case study
ANSYS/LS-DYNA based on ALE and Euler algorithm has the functions of heat transfer analysis, fluid and solid coupling analysis, static analysis and implicit analysis. It could quickly solve the dynamic nonlinear problems of high velocity impact in plane or space, explosion, and so on. ANSYS/LS-DYNA was adopted to simulate the impact of drift ice on tunnel. The drift ice is floating on the water surface. So it just need consider the horizontal loads such as wind drag force and flow drag force, while ignoring vertical loads such as gravity, buoyancy force, and so on. The parameters of tunnel and drift ice are shown in Table 1. The horizontal load is mainly expressed by drift ice initial velocity. When the drift ice with different size and thickness impacts the tunnel surface with different velocity, the impact force is different. The impact force can be simulated by modelling.
Physical property parameters of tunnel and drift ice
The automatic single contact (ASSC) was chosen as the contact type. Output control was conducted after assuming initial conditions of initial velocity, contact type and boundary condition, and then the K file was got. In the input calculation, the dynamics was used to analyze command stream file which is called K file in ANSYS, Fig. 3 shows part of K file. In the calculation, 3D SOLID164 solid element was adopted in the drift ice and tunnel model which are linear elastic material. The unit division adopted mesh mapping generation method as unit division. We consider the tunnel lining surface as the dominative surface and drift ice contact surface as the subsidiary surface. The parameters are shown in Table 1, where the parameters of ice were from the relation between the ice elastic modulus and ice in the references [11,12].
2.3 Functions and Effects of Drift Ice on Tunnel
The impact between the drift ice and tunnel is a complex collision problem between ice structures with many instability factors such as dynamic coupling, fluid-solid coupling, and so on. The main influence factors are velocity, plane size, thickness, shape, tensile strength, compressive strength, elastic modulus, flexibility of structure, friction, impact angle, and so on. This study mainly researches the impact force of drift ice on diversion tunnel surface in different velocity, different drift ice plane size, different drift ice thickness, and so on. In the calculation, other factors were ignored to reduce complexity. The tunnel finite element model is shown in Fig. 4
Tunnel finite element model.
Maximum impact force stress nephogram (v = 0.5 m/s).
2.3.1 Effect of drift ice velocity on impact force
The impact forces of drift ice crashing tunnel are different with different drift ice velocity. Set 2 m as water depth and 2 × 0.5 × 2 [TeX:] $$\mathrm{m}^{3}$$ as drift ice size, and calculate the drift ice impact force when the drift ice velocities are 0.5, 0.8, 1.0, 1.5, 1.8, 2.0, 2.3, 2.5, 2.8, 3.0 m/s, respectively. Within the damage permissible limits, the calculation results of the drift ice impact force are shown in Fig. 5.
Different velocities produce different stress and impact force when the drift ice crashing the tunnel, the flow velocity-maximum stress relationship graph and flow velocity-maximum impact force relationship graph are shown in Fig. 6.
(a) Flow velocity-maximum stress and (b) flow velocity-maximum impact force relational graph.
2.3.2 Effects of drift ice plane size on impact force
In order to research the influence of drift ice plane size on impact force, set 2 m as water depth, v = 5 m/s as drift ice velocity, 0.5 m as ice thickness. When the drift ice plane sizes are 0.5×0.5×2 m3, 1.0×0.5×2 m3, 1.5×0.5×2 [TeX:] $$\mathrm{m}^{3}$$, 2×0.5×2 [TeX:] $$\mathrm{m}^{3}$$, 2.5×0.5×2 m3, which are 1, 2, 3, 4, 5 times of 0.5×0.5 [TeX:] $$\mathrm{m}^{2}$$, respectively. Within the damage permissible limits, the calculation results of the drift ice impact force are shown in Fig. 7.
Different drift ice plane sizes produce different stresses and impact forces when drift ice crashing tunnel, the drift ice plane size-maximum stress relationship graph is shown in Fig. 8 and drift ice plane size-maximum impact force relationship is shown in Fig. 9.
Figs. 8 and 9 show that the collision stress relation curve is accordance with the collision impact force relation curve. When the drift ice has small plane size, the collision stress increases with the drift ice plane size, presenting linear relationship.
Maximum impact force stress nephogram in different drift ice plan size.
Drift ice plan size-maximum stress relational graph.
Drift ice plan size-maximum impact force relational graph.
2.3.3 Effects of drift ice thickness on impact force
The drift ice thicknesses change, the impact force of drift ice crashing tunnel changes accordingly. When other parameters are fixed, and the drift ice thicknesses are 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 m, respectively. Within the damage permissible limits, the calculation results of the drift ice impact force are shown in Fig. 10.
In different drift ice thicknesses, the relationship between maximum stress and thickness is shown in Fig. 11.
Fig. 12 shows that when other parameters are fixed, the impact force increases with the drift ice thickness when the thickness is less than 0.9m, and the stress value has a little change when the thickness is larger than 0.9m. When the drift ice thickness is small, the impact force increases with the thickness, presenting approximate linear relationship.
Maximum impact force local stress nephogram in different drift ice thickness.
Drift ice thickness-maximum stress relational graph.
Drift ice thickness-maximum impact force relational graph.
Above all, the tunnel lining surface would have different degrees of distortions when drift ice crashing tunnel lining in different drift ice velocity, different drift ice plane size, different drift ice thickness and so on. The distortions have little influence on the tunnel stability. But the tunnel lining surface breaks, and the tunnel lining surface would fall off by long time water erosion. It would destroy structure strength and stability. In the practical engineering, it would pay attention to this harm and adopt suitable protection measures.
2.4. Ice Mechanics Model Test
2.4.1 Test design
This study conducts indoor model test to research the impact mechanism of ice on tunnel in thawing period. In the indoor model test, the impact force of drift ice on diversion tunnel were tested on different drift ice plane size, different drift ice thickness, different drift ice velocity and so on. The model test state schematic diagram is shown in Fig. 13. Fig. 14 is the pump which used to pump water to other devices.
Model test state schematic diagram.
In this test, the main devices are a 2 m×1 m×1 m water tank which is used to store the water circulating in the system and a water pump which is the dominant body drives water from water tank to water channel, rises water level and operates circularly. They are shown in Fig. 15.
Test steps
The test has 15 steps:
Step 1: Fill the water tank with water.
Step 2: Open the main power and make water in the water tank rush into water channel by water pump.
Step 3: Put the ice blocks of 7 cm×3.5 cm×7 cm in water channel.
Step 4: Adjust flow velocity by rotary screw.
Step 5: The flow velocity is 0.5 m/s measured by flow meter.
Step 6: Measure the depth by water level probe and let the water depth be 7.2 cm.
Step 7: Stick 10 plexiglass strain gauges on the water line location of each side of tunnel model surface.
Step 8: Make the transient strain tester be close.
Step 9: Connect the strain gauges with transient strain tester by test wires and correspond each strain gauge with the joint in the transient strain tester.
Step 10: Turn on the transient strain tester and adjust to manual mode.
Step 11: Read the strain value which is the reading of the transient strain tester in sequence, along with the ice striking the tunnel surface.
Step 12: Record data, and process them.
Step 13: Change the conditions that set the flow velocity is 0.8, 1.0, 1.2, 1.5, 1.8, 2.0, 2.3, 2.5, 2.8, and 3.0 m/s, then repeat Step 4 to Step 12.
Step 14: Change the conditions that set the flow velocity is 1 m/s and set the ice blocks size to 7 cm×1 cm×7 cm, 7 cm×1.8 cm×7 cm, 7 cm×2.9 cm×7 cm, 7 cm×3.5 cm×7 cm, 7 cm×4.3cm×7 cm, 7 cm×5.5 cm×7 cm, respectively, then repeat Step 4 to Step 12.
Step 15: Change the conditions that set the flow velocity is 1 m/s and set the ice blocks size to 3.5 cm×1.8 cm×7 cm, 5.5 cm×1.8 cm×7 cm, 7 cm×1.8 cm×7 cm, 9 cm×1.8 cm×7 cm, respectively, then repeat Step 4 to Step 12.