** Damage Mechanism of Drift Ice Impact **

Li Gong* , Zhonghui Wang* , Yaxian Li* , Chunling Jin* and Jing Wang*

## Article Information

## Abstract

**Abstract:** The ice damage occurs frequently in cold and dry region of western China in winter ice period and spring thaw period. In the drift ice condition, it is easy to form different extrusion force or impact force to damage tunnel lining, causing project failure. The failure project could not arrive the original planning and construction goal, giving rise to the water allocation pressure which influences diversion irrigation and farming production in spring. This study conducts the theoretical study on contact-impact algorithm of drift ices crashing diversion tunnel based on the symmetric penalty function in finite element theory. ANSYS/LS-DYNA is adopted as the platform to establish tunnel model and drift ice model. LS-DYNA SOLVER is used as the solver and LS-PREPOST is used to do post-processing, analyzing the damage degrees of drift ices on tunnel. Constructing physical model in the experiment to verify and reveal the impact damage mechanism of drift ices on diversion tunnel. The software simulation results and the experiment results show that tunnel lining surface will form varying degree deformation and failure when drift ices crash tunnel lining on different velocity, different plan size and different thickness of drift ice. The researches also show that there are damages of drift ice impact force on tunnel lining in the thawing period in cold and dry region. By long time water scouring, the tunnel lining surfaces are broken and falling off which breaks the strength and stability of the structure.

**Keywords:** Damage Mechanism , Drift Ice , Impact

## 1. Introduction

Northwest China is located in the high, dry and cold region whose temperature is low and drift ice season is long in winter. Meanwhile, there is serious water shortage in this region, especially in part of Gansu Province, even the drinking water is in short supply. In order to ease the severe water shortage situation, the country successively invests and constructs a large numbers of water diversion projects for drinking and farm irrigation, such as the Datong River to Qinwangchuan region water diversion project, Tao River water diversion project, and so on. And the prepared west line of South-to-North [1] water diversion project is unprecedented. In winter operation period, the drift ice or ice berg would occur in the diversion projects under low temperature. They may cause serious impacts on water diversion project, including the influence of ice cover and ice berg on diversion capacity, the impact of ice expansion force on channel, the damage of drift ice on diversion tunnel, the impact of ice on channel operation and project security, and so on.

Domestic and foreign scholars have studied the evolution of ice conditions during the icing period, Domestic scholar Wang et al. [2] have conducted experimental research on the accumulation process of ice plugs by means of a sink test, Zhao et al. [3] conducted a study on the formation process of the ice and the formation of the ice in the Shensifenzi Bend, Wang et al. [4] introduced the influence of thermodynamic processes on the simulation of the sea ice and the elimination of the sea, and Shi et al. [5] based on the ice creep test, a nonlinear Burgers' sea ice model considering damage failure was established. The foreign scholar Jimenez et al. [6] established an updated-Lagrangian damage mechanics formulation for modeling the creeping flow and fracture of ice sheets, Marchenko and Cole [7] compared wave energy dissipation caused by three physical processes, Shortt and Sammonds [8] established the empirical scaling relations for sea ice mechanics by studying the micromechanics of ice. Jin et al. [9]proposed a gray GMP-Verhulst combination forecasting model identifies the change law of risk fluctuation. Gong et al. [10] used the symmetric penalty function in the finite element contact-impact algorithm to conduct the theoretical study on collision simulation problem between drift ice and water diversion tunnel.

In order to research the drift ice impact on the diversion tunnel in thawing period in long distance diversion project in cold and dry region, this study adopted theoretical analysis, numerical modelling and test research method. It intensively researched the impact force of the interaction between drift ice and diversion tunnel lining, using theoretical analysis, numerical simulation and testing research methods. The software of ANSYS/LS-DYNA was adopted as the platform to establish finite element model when the drift ice crashing diversion tunnel. The unilateralism slope experimental facility was also built to simulate evolution processes of drift ice crashing diversion tunnel. They reveal the impact influence regularity of drift ice on diversion tunnel on different working conditions. This study supplies theoretical support and technical guarantee for water conveyance project security in cold and dry region of western China.

## 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):

##### (3)

[TeX:] $$\ddot{u}_{j}=\frac{\partial^{2} u_{i}(X, t)}{\partial t^{2}}=\frac{\partial^{2}\left[x_{i}(X, t)-X_{i}\right]}{\partial t^{2}}=\frac{\partial^{2} v_{i}(X, t)}{\partial t^{2}}$$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).

##### (5)

[TeX:] $$D_{i j}=\frac{1}{2}\left(\frac{\partial v_{i}}{\partial x_{j}}+\frac{\partial v_{j}}{\partial x_{i}}\right)$$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

##### (7)

[TeX:] $$\int_{V} \delta v_{i}\left(\frac{\partial \sigma_{i j}}{\partial x_{j}}+\rho b_{i}-\rho \ddot{u}_{i}\right) d V=0$$where [TeX:] $$\delta v_{i}(X)$$ is the virtual velocity which is solved by Eq. (8).

##### (8)

[TeX:] $$\delta v_{i}(X) \in R_{0} \quad R_{0}=\left\{\delta v_{i}\left|\delta v_{i} \in C^{0}(X), \delta v_{i}\right|_{4}=0\right\}$$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:

##### (9)

[TeX:] $$\int_{V} \frac{\partial\left(\delta v_{i}\right)}{\partial x_{i}} \sigma_{j i} d V-\int_{V} \delta v_{i} \rho b_{i} d V-\int_{A_{i}} \delta v_{i} \bar{t}_{i} d A+\int_{V} \delta v_{i} \rho i i_{i} d V=0$$In this way, all particles velocity, acceleration, virtual velocity and strain rate could be expressed as the follow equation set:

##### (10)

[TeX:] $$\left\{\begin{array}{l} {\dot{u}_{i}(X, t)=N_{I}(X) \dot{u}_{i I}(t)} \\ {\ddot{u}_{i}(X, t)=N_{I}(X) \ddot{u}_{i I}(t)} \\ {D_{i j}=\frac{1}{2}\left(\frac{\partial \dot{u}_{i}}{\partial x_{j}}+\frac{\partial \dot{u}_{j}}{\partial x_{i}}\right)=\frac{1}{2}\left(\dot{u}_{i 1} \frac{\partial N_{I}}{\partial x_{j}}+\dot{u}_{j I} \frac{\partial N_{I}}{\partial x_{i}}\right)} \\ {\delta v_{i}(x)=N_{I}(x) \delta v_{i I}} \end{array}\right.$$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.

##### 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.

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

##### 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.

##### 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.

##### 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.

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.

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.

## 3. Results

##### 3.1 Effects of Flow Velocity on Impact Force

The drift ice velocity changes, the impact force and stress value on the tunnel surfaces changes accordingly, showing in Table 2.

By analyzing large number of test data, it is found that the test results are close to the finite element simulation results. Fig. 16(a) shows that the impact force variations are approximately the same in the software simulation results, test observation calculation results and the results of y = 10x. Therefore, the finite element calculation method is feasible. The impact force of drift ice increases with the flow velocity, appearing linear distribution.

##### 3.2 Effects of Drift Ice Size on the Impact Force

The drift ice sizes are different, the impact force and stress on tunnel surfaces are different. The values are shown in Table 3.

The comparison of test values and model simulation values are shown in Fig. 16(b). The test values are close to model simulation calculation values. The impact force increases with the drift ice size.

##### 3.3 Effects of Drift Ice Thickness on the Impact Force

By the tests, the drift ice thickness is different, the impact force and stress on the tunnel surfaces are different. The values are shown in Table 4.

The comparisons of test observation results and model calculation results are shown in Fig. 16(c). Fig. 16(c) shows that the impact force variations are approximately the same in the software simulation results, test observation calculation results and the results of y = 6.7x + 0.2. Therefore, the finite element calculation method is feasible. Within a certain range, the impact force of drift ice increases with the drift ice thickness, appearing linear distribution.

## 4. Conclusions

The comparison results of software simulation and test observation show that their trends are almost the same. The impact force of drift ice on tunnel lining increases with the drift ice thickness, appearing linear distribution, which verifies the reliability of the following simulation results.

The researches show that there are damages of drift ice impact force on tunnel lining in the thawing period in cold and dry region. By long time water scouring, the tunnel lining surfaces are broken and falling off which breaks the strength and stability of the structure. With the increasing of the size, thickness and the flow velocity of the drift ice, the impact force is different.

In different ice velocity, different size, different thickness of the drift ice, the tunnel lining surfaces would be destroyed and deformed when the drift ice crashing the lining. When the drift ice thickness is less than 0.9m, the impact force increases with the drift ice thickness, and the relationship between drift ice size and maximum stress appears linear relationship. When the drift ice thickness is larger than 0.9m, the impact force changes a little.

In the practical engineering, the two-phase movement between the ice and water is complex which exists multi-factor coupling effects among viscous force, drag force and other factors. The impaction between drift ice and tunnel lining is irregular. It needs further researches in the choice of the impact angle in the simulation.

## Biography

##### Li Gong

https://orcid.org/0000-0002-4824-5109He received B.S. degree in North China University of Water Resources and Electric Power in 2000, M.S. degree in Northwest A&F University in 2007, and a Ph.D. degree in Lanzhou Jiaotong University in 2014. Now he is a professor at the School of Civil Engineering, Lanzhou Jiaotong University, Lanzhou, China.

## Biography

## Biography

## Biography

## References

- 1 K. L. Yang, "Review and frontier scientific issues of hydraulic control for long distance water diversion,"
*Journal of Hydraulic Engineering*, vol. 47, no. 3, pp. 424-435, 2016.doi:[[[10.13243/j.cnki.slxb.20150824]]] - 2 J. Wang, B. Zhang, P. Chen, T. Liu, "Experimental study of ice jam accumulation during freezing period,"
*Journal of Hydraulic Engineering*, vol. 47, no. 5, pp. 693-699, 2016.doi:[[[10.13243/j.cnki.slxb.20151048]]] - 3 S. Zhao, C. Li, C. Li, X. Shi, S. Zhao, "Processes of river ice and ice-jam formation in Shensifenzi Bend of the Yellow River,"
*Journal of Hydraulic Engineering*, vol. 48, no. 3, pp. 351-358, 2017.doi:[[[10.13243/j.cnki.slxb.20160721]]] - 4 K. Wang, P. Liu, S. Jin, N. Wang, Z. Yu, "Sea-ice growth and decay model of Bohai Sea based on thermodynamic process,"
*Advances in Water Science*, vol. 28, no. 1, pp. 116-123, 2017.doi:[[[10.14042/j.cnki.32.1309.2017.01.013]]] - 5 C. Shi, Y. Luo, Z. Hu, "Non-linear Burgers’ sea-ice model considering damage effects and its numerical application,"
*Engineering Mechanics*, vol. 35, no. 7, pp. 249-256, 2018.custom:[[[-]]] - 6 S. Jimenez, R. Duddu, J. Bassis, "An updated-Lagrangian damage mechanics formulation for modeling the creeping flow and fracture of ice sheets,"
*Computer Methods in Applied Mechanics and Engineering*, vol. 313, pp. 406-432, 2017.doi:[[[10.1016/j.cma.2016.09.034]]] - 7 A. Marchenko, D. Cole, "Three physical mechanisms of wave energy dissipation in solid ice," in
*Proceedings of the 24th International Conference on Port and Ocean Engineering under Arctic Conditions*, Busan, Korea, 2017;custom:[[[-]]] - 8 M. W. Shortt, P. R. Sammonds, "Experiments on the micromechanics of ice using scanning electron microscopy," in
*Proceedings of the 25th International Conference on Port and Ocean Engineering under Arctic Conditions*, Delft, Netherlands, 2019;custom:[[[-]]] - 9 C. Jin, M. Wu, L. Gong, "Risk prediction of ice-jam disaster in Ningxia-Inner Mongolia reaches of the Yellow River based on grey Markov-GMP-Verhulst model,"
*Journal of Natural Disasters*, vol. 28, no. 2, pp. 82-91, 2019.custom:[[[-]]] - 10 L. Gong, Y. Li, C. Jin, "Numerical simulation and verification on impact damage mechanical property of drift ice on diversion tunnel,"
*Transactions of the Chinese Society of Agricultural Engineering*, vol. 34, no. 13, pp. 144-151, 2018.doi:[[[10.11975/j.issn.1002-6819.2018.13.017]]] - 11 T. L. Yu, Z. G. Yuan, M. L. Huang, "Experimental study on mechanical behavior of river ice,"
*Journal of Liaoning Technical University (Natural Science)*, vol. 28, no. 6, pp. 937-940, 2009.custom:[[[-]]] - 12 F. Zhang, G. Li, "Simulation of flow ice impact on dam body based on LS-DYNA,"
*Water Conservancy Construction and Management. 0.9m0.9m*, vol. 33, no. 2, pp. 19-21, 2013.custom:[[[-]]]