Molecular Dynamics Simulation and Analysis of Crystallization System in PVC Tunnel Drain-Pipe at Different Temperatures

Corresponding Author: Xuefu Zhang School of Civil Engineering, Chongqing Jiaotong University, Chongqing, 400074, China Email: zhangxuefu400074@126.com Abstract: Aiming at the problem that carbonate and calcium ions are easy to crystallize and precipitate in the drain-pipes of limestone tunnels and then adsorb on PVC pipes in Western China, a "solid-liquid interface model" for the interaction between PVC materials and ions in aqueous solution is established by using the molecular dynamics simulation software materials studio and the microscopic mechanism of crystallization system in PVC drain-pipes at different temperatures is studied. The results show that the selfdiffusion coefficient of water molecules increases with the increase of temperature. In the temperature range of 283-308 k, the self-diffusion coefficient of calcium ion and carbonate decreases first, then increases and decreases again with the increase of temperature. When T = 303 K, the selfdiffusion coefficient of both ions reaches the maximum value, which makes it easier to crystallize. With the increase of temperature, the binding energy between ionic solution system and PVC increases at first and then decreases. The research results and methods have important guiding significance for the prevention and control of crystal blockage of limestone tunnel drainage pipe.


Introduction
With the development of traffic planning in central and western regions of China, the importance of tunnel construction has been paid more and more attention. There are many limestone tunnels in Southwest China. The water in this limestone geological environment contains nearly saturated bicarbonate and calcium ions. The scaling ions flow into the drainage pipeline along with the underground water and form calcium carbonate crystal in the tunnel drainage system and deposit in the drain-pipe (Chen et al., 2019;Wu et al., 2019;Al Nasser and Al Salhi, 2015). It is easy to cause the blockage of tunnel drainage system, causing the rupture and leakage of tunnel lining, the argillization of surrounding rock structural plane weak interlayer due to water erosion and the cracking under freeze-thaw cycle of drainage pipeline and lining in alpine region.
Most of the research on tunnel drainage is mainly reflected in discovery and treatment (Lee et al., 2012;Jung et al., 2013) and the research on pipeline crystallization scaling is less. As early as the end of 1950s, Kern put forward the mass balance model, believing that scaling can be divided into deposition and removal and the difference between them is the amount of scaling (Kern, 1959). Then, (Taborek et al. 1972) established a prediction model based on the mass balance model and considered that the surface reaction was the main factor to control the scaling deposition. Hasson et al. (1968) studied the mechanism of CaCO3 scale deposition on the surface of heat exchanger by water turbulence containing dissolved scale component under nonboiling condition and proposed a CaCO3 scale deposition rate model. Since the 21st century, the research on crystallization of calcium carbonate in aqueous solution and pipeline scaling has entered the stage of numerical model research through geochemical calculation software such as PHREEQC. For example, (Tao et al., 2007) used PHREEQC to simulate and analyze the uranium leaching agent of groundwater in a high salinity area in China and determined the critical value related to hydrogeochemical conditions of crystallization precipitation through the basic conditions of crystallization, such as saturation index. Zhu analyzed the effects of pH, Ca 2+ and bicarbonate on the dissolution and precipitation of calcium carbonate. Chen studied the influence of activity product, solubility product and supersaturation index on the crystallization process through related crystallization experiments and PHREEQC. Some scholars also use molecular dynamics to analyze the adsorption behavior of molecular ions. Hu respectively studied the adsorption behavior of protein on the surface of PEG and PDMS polymer antifouling films. Huang et al. (2005) simulated and analyzed the interaction between Polyvinylidene Fluoride (PVDF), Polytrifluorochloroethylene (PCTFE), Fluororesin (F2314), Fluororubber (F2311) and 1,3,5-Triamino and 2,4,6-Trinitrobenzene (TATB) and obtained the binding energies between four fluoropolymers and TATB. The binding energies of four fluoropolymers with TATB were obtained. Tian (2008) discussed the scaling process of Fe (001) crystal face and (Kong and Wang, 2016) analyzed the adsorption behavior of Cu (II) on hydroxylated kaolin (001) surface in aqueous environment by molecular dynamics.
It can be seen that the microscopic study of tunnel drain-pipe wall on the crystallization of scaling ions in aqueous solution on the pipe wall surface is still lacking, especially on the interaction between PVC materials and ions in aqueous solution. In general, a large number of systematic experiments are needed to study. The experiments usually need complete equipment, timeconsuming and laborious. Moreover, it is not possible to study the related properties between the two interfaces from the microscopic point of view (Kirkham et al., 2020;Chihi et al., 2020). Therefore, molecular dynamics method is used in this study to study the adsorption properties of PVC materials to ions in aqueous solution and the "solid-liquid interface model" is used to explain or predict the formation of crystallization system in PVC drainage pipe. Meanwhile, the micro mechanism of crystallization in limestone tunnel drain-pipe is studied. The results and methods have important guiding significance for the prevention and control of crystallization blockage of limestone tunnel drain-pipe.

Construction of Interface Model between PVC Layer and Solution
The process map of construction and method of numerical model are shown in Fig. 1. For the construction of PVC layer polymer, the larger the polymerization degree of molecular chain is, the closer to the real model and the more accurate the simulation results are, but this requires higher computer configuration and the longer the actual simulation process will take (Emery et al., 2020;Lewington et al., 2020;Schmidt et al., 2020). Therefore, the appropriate degree of polymerization is needed in the actual model calculation. In this study, considering the accuracy and time of calculation, the polymerization degree of 10 is selected. PVC chain is established by materials visualizer module in Materials Studio software. The molecular formula is shown in Fig. 2a and the single chain model is shown in Fig. 2b. Through Modules->Amorphrouscell->construction module, the amorphous three-dimensional structure of 10 chains is established. The size of the three-dimensional crystal cell structure is 19.5172×19.5172×19.5172 Å 3 and the crystal cell structure is shown in Fig. 3a. Similarly, the visualizer and amorphous cell modules in Materials Studio software are used to construct calcium carbonate aqueous solution. The calcium carbonate aqueous solution system constructed in this study contains 300 water molecules, 3 carbonate ions and 3 calcium ions. The volume of the system is 19.5172×19.5172×24.8605 Å 3 . To compare the interaction between calcium ion and carbonate ion with polyvinyl chloride layer, the two ions are placed at the bottom of the solution near PVC. After the completion of the construction, it is also optimized by molecular mechanics and the final interface model of water solution and PVC layer is shown in Fig. 3b.

Determination of Simulation Parameters
The discovery module suitable for Molecular Dynamics (MD) simulation in Materials Studio software is used. Using high precision universal compass force field, all atom coordinates of PVC layer are fixed. Since the pressure of the system is not the key factor, the canonical ensemble (NVT) is adopted. Firstly, the NVT ensemble and velocity scale temperature control method are used to simulate the dynamics of the system at 100 ps to make the system reach the equilibrium state. Secondly, the kinetic simulation is carried out in the NVT ensemble and Andersen constant temperature hot bath. The time step of the process is 1 fs and the simulation time is 200 ps. The trajectory of the system is recorded once every 10000 steps and the simulated temperature is 283-308 k. The simulation is conducted once every 5K and the truncation radius is 9.5 Å.

Determination of System Equilibrium
The system equilibrium is judged by temperature and energy balance. The accuracy of the simulation is characterized by observing the energy convergence parameter, the ratio of total energy fluctuation value and kinetic energy fluctuation value (R). When the energy convergence parameter ≤0.001, R≤0.001, the calculation results are reliable. The equilibrium process of the system less than 298 K is illustrated as an example. Figure 4a shows the energy output curve of the equilibrium process and Fig. 4b shows the temperature output curve of the equilibrium process. It can be seen from the energy curve that the change of potential energy and non-bond energy of the system with time has been flattened, indicating that the energy of the system has reached the equilibrium state. It can be seen from the temperature curve that the temperature fluctuates 10% above and below 298 K, indicating that the system temperature has reached the equilibrium state and the determination of equilibrium at other temperatures is also based on this method.

Analysis of Self-Diffusion Coefficient
To further study the dynamic characteristics of the system, the self-diffusion coefficients of the system at different temperatures are analyzed and calculated. The self-diffusion coefficient D is the physical quantity measuring the motion of a single colloidal particle surrounded by the same uniformly distributed particles (Tian, 2008). There are many methods to calculate the self-diffusion coefficient. In this study, the selfdiffusion coefficient of the particle is calculated according to the mean square displacement of the particle. Based on the mean square displacement expression and Einstein's diffusion law, the value of diffusion coefficient D can be obtained. Figure 5 shows the mean square displacement curves of water molecules at different temperatures obtained by molecular dynamics simulation at 283k-308k. When it is less than 10ps, it is a transition period from an unstable state to a stable state and the mean square displacement curve after 10ps is an approximately straight line. It can be found that with the increase of temperature, the slope of the mean square displacement curve increases, thus the self-diffusion coefficient of water molecules also increases. This is because the kinetic energy of water molecules increases with the increase of temperature and the hydrogen bonds between some water molecules are destroyed, resulting in an increase in the number of free single water molecules. Thus the self-diffusion coefficient of water molecules is increased, that is, the activity of water molecules is enhanced. At the same time, with the increase of the activity of water molecules, the contact probability of calcium ions and carbonate ions will be reduced. Therefore, the increase of temperature will inhibit the crystal growth to a certain extent. In the same way, the self-diffusion coefficients of calcium ion and carbonate ion can be obtained as shown in Fig. 6. In the temperature range of 283-308k, the selfdiffusion coefficients of calcium ion and carbonate first decrease, then increase and decrease again. The selfdiffusion coefficients of calcium ion and carbonate ion change obviously at different temperatures, which indicates that temperature has a great influence on the self-diffusion coefficient of calcium ion and carbonate ion. When t = 303 k, the self-diffusion coefficients of calcium ion and carbonate ion reach the maximum value. At this time, the activity of the two ions is large and they are easy to react to form ion pairs and crystals under certain conditions. Moreover, it can be seen that the temperature is an inflection point. If it is higher than this temperature, the diffusion coefficient tends to decrease. The simulation results of the self-diffusion coefficient in the temperature of 283-308 k are similar to those obtained in other references (Wang, 2010;Wang et al., 2010), which verifies the correctness of the results.

Analysis of Interaction between Systems
The interaction between PVC and scaling ion solution system is analyzed by simulation. If the interaction is very strong, the scaling ion solution system is easy to adhere to the PVC layer and the pipe wall is easy to crystallize. To calculate the adsorption energy of PVC layer, the interaction energy of the system is defined as ΔE and the binding energy Ebinding as the opposite number of interaction energy ΔE. Details as follows: In Eq. (1), Etotal is the total energy of system, Ep is the single point energy of PVC layer and Eo is the single point energy of calcium ion, carbonate ion and water molecule after interaction (Tian, 2008). The calculated energy values of the system under different temperatures are shown in Table 1. The change trend of binding energy with temperature is shown in Fig. 7.
It can be seen from Table 1 that the interaction energy is negative, which means that the adsorption of scaling ion solution system on PVC surface is a spontaneous process and a relatively stable system can be formed. Compared with (Tian, 2008), the adsorption energy of PVC for calcium ion and carbonate ion system is far less than that of calcite crystal surface. Figure 6 is a graph of the binding energy between the ionic solution system and PVC changing with temperature. Combined with Table 1 and Fig. 6, it can be seen that within the range of simulated temperature, the binding energy between scale ion solution system and PVC increases at first and then decreases with the increase of temperature. When the temperature T = 283 K, the binding energy is the minimum and Ca 2+ and CO3 2are the most difficult to deposit and crystallize on the surface of PVC pipe wall. When temperature T = 298 K, the binding energy between ionic solution system and PVC reaches the maximum and Ca 2+ and CO3 2are easy to deposit on the surface of PVC pipe wall.

Conclusion
A two-layer model system of PVC layer and aqueous solution containing calcium ion and carbonate ion was established and molecular dynamics simulation of the system at different temperatures was carried out. The interaction energy, self-diffusion coefficient and interaction relationship between PVC layer and aqueous solution system were analyzed. It can be concluded that: (1) With the increase of temperature, the self-diffusion coefficient of water molecules increases. The increased activity of water molecules will reduce the contact probability of calcium ions and carbonate ions. Therefore, the increase of temperature will inhibit the crystal growth to a certain extent (2) In the temperature range of 283-308k, the selfdiffusion coefficients of calcium ion and carbonate first decrease, then increase and decrease again. When T = 303K, the self-diffusion coefficient of calcium ion and carbonate ion reaches the maximum value. At this time, the activity of the two ions is relatively large and it is easy to react to form ion pairs and crystals under certain conditions (3) The adsorption energy of PVC for calcium ion and carbonate ion solution system is much less than that of calcite crystal surface. In a certain temperature range, the binding energy between the ionic solution system and PVC increased first and then decreased with the increase of temperature When the temperature T = 283k, the binding energy is the minimum and Ca 2+ and CO3 2are the most difficult to deposit and crsystallize on the surface of PVC pipe wall. When temperature T = 298K, the binding energy between ionic solution system and PVC reaches the maximum and Ca 2+ and CO3 2are easy to deposit on the surface of PVC pipe wall.
Due to the consideration of the accuracy and calculation time of the model results in this study, the length of the degree of polymerization and the molecular concentration of calcium carbonate solution are limited. In the future, we will conduct further in-depth research on this problem, so that molecular dynamics simulation in this field can more truly approximate the field situation.