A non-linear flow model for the flow behavior of water inrush induced by the karst collapse column

Water inrush induced by the karst collapse column (KCC) is a great threat to coal mine safety. In this study, a non-linear flow model that couples three flow types is built based on flow transition from laminar flow in the aquifer to turbulent flow in the mine roadways during the process of water inrush induced by KCC. The proposed model couples Darcy flow, Forchheimer flow, and turbulent flow, and is then used to simulate the flow behavior of water inrush induced by KCC. In particular, the “3.1” water inrush incident from the coal seam floor in the Luotuoshan coal mine, China, is numerically investigated. The numerical results show that with the increase of the inrush flow rate, Forchheimer flow in the water inrush channel is first controlled by viscous resistance, then affected by both viscous resistance and inertial resistance, and finally controlled by inertial resistance. Therefore, water inrush induced by KCC is a dynamic process with a transition from laminar to turbulent. The Forchheimer equation proved to be applicable in describing the high-velocity non-linear flow, and can also reflect the intermediate state of the flow translation from laminar flow in the aquifer to turbulent flow in the roadway during the water inrush process.


Introduction
The karst collapse column (KCC) is one special type of karst collapse in the permo-carboniferous coal eld of North China. The KCC is composed of a fractured rock mass of coal measure strata. Water inrush induced by KCC can be sudden and catastrophic, 1,2 which threatens the safety of coal mines. KCCs are widely dispersed in almost 45 coal mines attached to 20 coal elds in Shanxi, Shandong, Henan, Jiangsu, Hebei, and Shaanxi provinces of China. 2,3 Many Chinese coal seams lie above the carbonate karst rocks. The total number of Chinese coal mines has increased to over 3000. 4 Long tunnel construction projects are also threatened by karst water inrush. According to the relevant statistical data, large-scale karst water inrushes have taken place in almost 50% of the 17 karst tunnels located in southwest and middle-south China since 1989. 5 Water inrush induced by KCC indicates that the rock mass in the KCC has hydraulically connected a karst aquifer and the mine roadway. The relationship between the velocity and pressure gradient of the inrush in a fractured rock mass generally does not satisfy the Darcy equation, and is distinctly nonlinear. [6][7][8][9][10] It is difficult to describe the non-linear ow pattern based on the Darcy equation, which is applicable in linear ow only. Therefore, it is of great theoretical and practical signicance to build a nonlinear ow model for the identication of ow mechanisms and the reasonable forecast of water inrush.
Two main equations are used to describe the non-linear relationship between the pressure and ow rate of uids in porous media: the Izbash equation and the Forchheimer equation (see f.e.). 11,12 The Forchheimer equation was rst proposed based on experiment and then later demonstrated by theoretical inference. 13,14 The rst-degree term of the Forchheimer equation is associated with viscous resistance, and the quadratic term is associated with inertial resistance. When the ow rate is small enough, the viscous resistance is the main factor affecting the ow behavior. However, when the ow rate is high, inertial resistance is the main factor. Compared with the Izbash equation, the Forchheimer equation has clear theoretical and physical background, which can describe the high-velocity non-linear ow in porous media (with large porosity) and in fractured media. [15][16][17][18] Water inrush induced by KCC in mining engineering shows the characteristics of Forchheimer ow. Miao et al. 19 and Ma et al. 6 conducted seepage experiments on broken rocks and found that the seepage in the broken rock obeys the Forchheimer equation. Sedghi-Asl et al. 20 conducted seepage experiments on different aggregate size gravels and found that the non-Darcy factor decreases as the aggregate size increases based on the Forchheimer equation. Different types of porous media were also examined to experimentally study the parameters of the Forchheimer equation and these results were systematically summarized by Moutsopoulos et al. 21 In recent years, numerical simulation has been widely used to study high-velocity non-Darcy ow behavior in sand and gravel riverbeds, 22 rocklls, 23 and earth and rockll dams. 24 Basak et al. 25 presented the effect of non-linearity in the ow response on the discharge characteristics and pressure distribution of a non-penetrating well in a semi-innite medium incorporating the Forchheimer equation. Wang et al. 26 also performed similar studies. These models are all based on the Forchheimer equation. Moreover, Xu et al. 27,28 proposed an innovative simulation method to study the ow state evolution laws in the karst regions by coupling Darcy's Law, Brinkman equations and the incompressible Navier-Stokes equations.
Groundwater travels from an aquifer though a KCC into a mine roadway during the process of water inrush induced by the KCC. The history of the pressure and velocity of each ow eld is a time-varying physical process and hence, the three ow elds are inseparable both in time and space. 29 Although laboratory tests 20,22 and theoretical analysis 13,22,30 are precise and reliable in revealing the non-linear ow mechanisms of Forchheimer ow in fractured rock masses, the ow mechanisms of Forchheimer ow for large-scale engineering problems still cannot be quantitatively analyzed, particularly for the history of ow and pressure of different ow elds during the inrush process. Therefore, it is important to build a coupled non-linear ow model considering the three ow elds: Darcy ow in the aquifer, high-velocity non-linear ow in the KCC and turbulent ow in the roadway during the entire ow process of water inrush.
In this study, in order to explain nonlinear ow behavior during the entire process of water inrush with better accuracy, the essence of ow transition was rst veried by experiments. Then, a non-linear ow model coupling three ow elds was established based on the ow transition. The proposed model was nally used for numerical simulation of the process of water inrush in the Luotuoshan coal mine, China. In addition, the effect of the KCC permeability on ow behavior was discussed and some main conclusions were drawn.

Experiment on the flow translation of water inrush
Engineering practice shows that a water inrush disaster in a coal mine is the physical process of groundwater bursting from the aquifer into the mine roadway, which is induced by the combination of mining stress and high water pressure. 1,19,31 The quantitative change of the inow shows ow transitions from laminar ow in the aquifer to turbulent ow in the roadway. Zhao et al. 32 postulated that water inrush from a bearing cave is induced by the failure of a rock pillar. Aer the failure of the rock pillar, water bursts out from the conned karst cave to form pipe ow. In physics, Darcy laminar ow is viscous resistance-based, and Navier-Stokes turbulent ow is inertial resistance-based. The state of Forchheimer ow is relatively complex and is based on both viscous resistance and inertial resistance. The classical theory of uid mechanics uses the relationship between the Reynolds number and the Fanning friction coefficient to judge the transition between the above three states. 33 In this section, a seepage experiment is conducted to study the ow transition and to provide an experimental basis for the establishment of the coupled non-linear ow model.
The testing system of non-linear ow in porous media, independently developed by Northeastern University, China, was employed to study the non-linear ow behavior of different porous media, as shown in Fig. 1a. The sample, 320 mm long and 60 mm in diameter, was saturated aer being packed with sand. A constant water pressure was maintained during the experiment. The water entered from the bottom of the sample, and exited from the top of the sample. Piezometers were used to determine the pressure change along the column during the experiment and to calculate permeability. An electronic balance was used to measure water inow. Data collection was carried out automatically by a computer with a 3 s sampling interval. The pre-set range of the hydraulic gradient was 3-36.  34 Due to the inhomogeneity of the sands, the mean particle diameter was used as the characteristic length to calculate the Reynolds number and the Fanning friction coefficient. 33 A relationship between the Reynolds number, Re, and the Fanning friction coefficient, f, was obtained and shown in Fig. 2. It can be observed that the ow in ne particle size sands satises Darcy's law. However, the relationship between Re and f gradually deviates from linearity with an increase in particle size, which indicates that the relationship deviates from Darcy's law and gradually transforms to Darcy-Forchheimer ow. When Re is signicantly high, the Fanning friction coefficient does not decrease as the Reynolds number increases, but increases gradually. This indicates that the ow state is transitioning from laminar ow to turbulent ow, which has been conrmed by Tzelepis et al. 35 In summary, under constant water pressure gradient, the ow state in porous media with different particle size sands is not necessarily the same. There are two ow transitions during the entire process of the experiments: (1) the transition from linear laminar ow to inertial ow; (2) the transition from inertial ow to turbulent ow. The larger the particle size, the more easily non-Darcy ow occurs. Therefore, for water inrush induced by KCC, the ow in the KCC may be in the intermediate state between laminar ow in the aquifer and turbulent ow in the roadway, and inevitably experiences ow transition from laminar ow to turbulent ow. This has also been veried through some studies by Yang et al., 34 Ma et al., 6,7 and Zhang et al. 36 A ow model based on one ow state cannot reect the ow transition of water inrush.

Laminar ow in aquifer
Groundwater ow rate in an aquifer is usually small and the ow resistance is primarily viscosity dominated. The ow state of normal groundwater is laminar ow. The uid velocity and pressure gradient satisfy the linear Darcy's Law, 33 which can be expressed as where v D is the ow velocity (m s À1 ), p D is the pressure (Pa), k D is the permeability (m 2 ), and m is the dynamic viscosity (Pa s). The continuity equation 33 is where r is the density of groundwater (kg m À3 ), f D is the porosity of aquifer, Q D is the source (sink) term (kg m À3 s À1 ), and t is the time (s).

Forchheimer ow in KCC
The Forchheimer equation is widely used in the study of ow characteristics of earth and rockll dams. Flow experiments conducted on fractured rocks indicate that the ow exhibits evident non-linear behavior and the relationship between the velocity and the pressure gradient effectively obeys the Forchheimer equation. Flow in fractured rocks is still laminar at a certain ow rate. The non-linear ow momentum equation of incompressible uid in fractured rocks based on the Forchheimer equation 37 can be written as The continuity equation 33 is given as The relation between permeability and the non-Darcy factor can be given as 38,39 where c a is the acceleration coefficient, v F is the ow velocity (m s À1 ), p F is the pressure (Pa), k F is the permeability (m 2 ), b F is the non-Darcy factor (m À1 ), c is the Forchheimer coefficient, which can be obtained from seepage experiments, Q F is the source (sink) term (kg m À3 s À1 ), and f F is the porosity of KCC.

Turbulent ow in roadway
The groundwater ow state in a mine roadway is turbulent ow, which satises the viscous Newtonian uid Navier-Stokes equation, 40,41 and can be expressed as 8 < : where v NS is the ow velocity (m s À1 ), p NS is the pressure (Pa), and F is the body force (N m À3 ).

Boundary conditions of the adjacent ow
During the entire process of groundwater travelling from the aquifer to the KCC and discharging into the roadway, both the ow rate and the pressure are continuous. Therefore, the transition boundary condition can be dened as follows: On the boundary between the aquifer and the KCC & On the boundary between the KCC and the roadway & The above non-linear ow model is the combination of the ow in the aquifer, KCC, and roadway. This model can reect the ow transitions from laminar ow in an aquifer to turbulent ow in a roadway.

Numerical methods and verication
The numerical solution of Darcy and Navier-Stokes is easy to implement and can be directly called in COMSOL. The nonlinear Forchheimer equation is solved by applying the Galerkin nite element method and simple linearization method. 42 A more detailed analysis of numerical astringency while solving the equations has been presented in the literature. 27 The nite element calculation program is compiled in Matlab, and the numerical solution of the coupled ow eld can be realized through iterative solution using the module of COMSOL with Matlab.
A constant pressure is set for the inlet boundary of the aquifer, and a velocity, equaling the ow velocity of KCC at the interface, is set for the outlet boundary of the aquifer. As for the KCC, a pressure is set for the inlet boundary, which equals the pressure of the aquifer at the interface, and a velocity is set for the outlet boundary, which equals the ow velocity of the roadway at the interface. Furthermore, a pressure is set for the inlet boundary of the roadway, which equals the pressure of the KCC at the interface, and a constant pressure is set for the outlet boundary of the roadway. Specic details about the boundary conditions of the adjacent ow are provided in Table 1.
As the built-in modules of COMSOL, the Darcy equation and the Navier-Stokes equation have been widely used. Therefore, only the nite element method for the Forchheimer equation needs to be veried. Fig. 3 shows the results of the seepage experiment of porous media with a particle size of 2.00-2.36 mm, from which the permeability, k F , and the non-Darcy factor, b F , of the Forchheimer equation can be obtained, i.e., k F ¼ 6.465 Â 10 À10 m 2 and b F ¼ 1.735 Â 10 5 m 2 . Based on the two medium parameters, the corresponding relationship between the pressure gradient and the ow velocity is obtained through numerical simulation. Fig. 3 shows that the numerical results are in good agreement with the experimental tting results. Aquifer This journal is © The Royal Society of Chemistry 2018

A case study
As a case study of water inrush from a coal seam oor, a numerical simulation based on the proposed non-linear ow model was conducted to investigate the ow behavior of the "3.1" water inrush incident in the Luotuoshan coal mine, China.
The water inrush occurred in the no. 16 coal seam on the +870 level of the Luotuoshan coal mine of the Shenhua Group on March 1, 2010. The maximum water inrush was up to 72 000 m 3 h À1 . In the incident, the mine was ooded, 32 miners died, and 7 miners were injured. The direct economic loss was about RMB 48 million. A hydrogeological investigation and a pumping test indicated that the water inrush source was the Ordovician limestone karst aquifer with abundant water content. The average thickness of the aquifer is about 23 m, and the Ordovician aquifer water pressure is about 4.1 MPa. Hydrogeological investigation and similar material model tests indicate that the water-conducting channel is a developing KCC. The speculated section of the KCC is shown in Fig. 4. The area of the excavated cross section of the no. 16 coal seam is 19 m 2 . The water inrush occurred at the working face of the coal seam roadway.

Numerical model
According to the geologic and geometric structure of the water inrush induced by the KCC in the Luotuoshan coal mine, a two dimensional numerical model is established as shown in Fig. 5. The model comprises three parts: the limestone aquifer with an area of 100 m Â 23 m, the water-conducting channel of the KCC with a maximum width of 56 m, height of 21.5 m, and the roadway with an area of 56.2 m Â 3.8 m. According to the hydrological conditions, the boundary conditions of the model were set as follows. A constant hydraulic pressure, p ¼ 4.1 MPa, was set for the le and right boundaries of the aquifer. Constant atmospheric pressure, p ¼ 0.1 MPa, was set for the right exit of the roadway. At the interface of the aquifer and the KCC, p D ¼ p F and v D ¼ v F were set. At the interface of the KCC and the roadway, p F ¼ p NS and v F ¼ v NS were used. A boundary condition with zero ux was applied at the rest of the external boundary of the model. An initial water pressure, p ¼ 4.1 MPa, was set for the aquifer and the KCC. According to the eld test, the hydrodynamic parameters are shown in Table 2.
The accuracy of the nite element calculation is largely related to the quality of the nite element mesh. In general, the smaller the mesh size, the more accurate the calculation result and the more time-consuming the calculation. Thus, an error analysis was rst conducted to determine the most suitable mesh size in order to eliminate the effect of the nite element mesh on the calculation results as much as possible. Fig. 6 shows the results of groundwater inow and relative error under different maximum mesh size using k F ¼ 9.6 Â 10 À9 m 2 . As shown in Fig. 6, the calculation error increases with the maximum mesh size. Considering the calculation efficiency and the calculation error, the model is discretised into a mesh that contains 19 370 six-node triangle elements with a maximum mesh size of 0.6 m and a relative error of 1%. A time step of 2 s was set for the models and the total time was 180 s.

Numerical results
Numerical simulation of the water inrush was conducted using k F ¼ 9.6 Â 10 À9 m 2 . The results are shown in Fig. 7-10. Fig. 7 and 8 represent the distributions of the water velocity and pressure. Fig. 9 shows the distribution of the streamlines and the velocity vectors. The water inow history of the measuring point M is shown in Fig. 10. As shown in Fig. 7 and 8, the ow velocity and pressure changes continuously during the entire process. The ow state of the inrush in the KCC is Forchheimer ow, between Darcy laminar ow dominated by the viscous resistance and turbulent ow dominated by the inertial resistance. The ow velocity increases with time. Forchheimer ow is rst dominated by viscous resistance. Gradually, it is dominated by both inertial resistance and viscous resistance. Finally, Forchheimer ow is completely dominated by the inertial resistance. Therefore, the water inrush induced by a KCC is a dynamic process of ow changing from gradual to sudden. As  shown in Fig. 10, water ow with high pressure in the KCC rushes into the roadway continuously as the uid pressure in the KCC is translated into momentum of the uid. Therefore, the pressure on the entrance of the KCC drops gradually; conversely, the water inow increases gradually. As the conned aquifer boundary pressures remain constant, the pressure distribution in the KCC and aquifer will eventually attain equilibrium according to the pressure balance principle. In addition, the computed water ow at the outlet of the roadway will reach a steady state.

Sensitivity analysis of the KCC permeability
During the water inrush process induced by KCC, the KCC permeability constantly increases, due to the loss of sediment and broken stone particles under the effects of high pressure and high rate ow scouring. 7,19 Under certain water pressure conditions, the KCC permeability and the non-Darcy factor not only directly affect the relative weights of the viscous resistance and inertial resistance, but also decide the water pressure distribution in the KCC and the water inow in the roadway. There is a one-to-one correspondence between the permeability and non-Darcy factor. 21,22 Therefore, sensitivity analysis of the KCC permeability has implications for the elucidation of the mechanism of water inrush. 4.3.1. The method of discriminating the ow state. The Forchheimer number is the ratio of the quadratic term of the Forchheimer equation (the inertial term) and the rst-degree term (viscosity term), which could reect the non-linear degree of the non-linear ow. It was applied in discriminating non-linear ow by Ruth and Ma 43 and Zeng and Grigg: 44 Eqn (9) shows that when the ow velocity is sufficiently small (F o ( 1), the ow inertial resistance can be ignored compared with the viscous resistance. The ow state tends to be Darcy laminar, which shows weak nonlinearity. When the ow velocity is high enough (F o [ 1), the ow inertial resistance cannot be ignored compared with the viscous resistance. The ow state  6 The effect of mesh size on groundwater inflow and relative error. This journal is © The Royal Society of Chemistry 2018 tends to be turbulent, which shows strong nonlinearity. When the velocity is in a certain range (F o z 1), the inertial resistance is about the same as the viscous resistance. Neither of the resistances can be ignored, and the ow state is non-linear laminar controlled by both the inertial resistance and viscous resistance. 4.3.2. The effect of permeability on the ow state. In order to study the effect of the KCC permeability on the ow state, a variable n was dened as the ratio of the KCC permeability to the aquifer permeability, which can be expressed as n ¼ k F /k D . There is a one-to-one correspondence between each n and the permeability; the greater the value of n, the higher the KCC permeability. In this study, the variable n is set in a wide range (10 0 -10 3 ). For each magnitude of n, the water pressure and velocity on the ow path when the ow is stable is extracted and illustrated in Fig. 11. It can be observed that the uid pressure and velocity in the KCC are very sensitive to the KCC permeability. The permeability affects the non-Darcy effect and thus, affects the non-linear inertial resistance. Under the constant aquifer pressure boundaries, the higher the KCC permeability, the greater is the pressure drop in the aquifer. The ow velocity at the same location in the KCC increases as the KCC permeability increases. Because the KCC is shaped like a plug (the width of the bottom is obviously bigger than the top), according to the mass conservation of uids, the narrower the KCC channel, the greater is the ow velocity. Taking n ¼ 10 3 as an example, the highest velocity in the channel can be up to 0.13 m s À1 , and is located in the narrowest part of the channel. Fig. 12 shows the relationship between n and the non-Darcy factor, Forchheimer number, water inow and water pressures at the entrance of the KCC (measuring point M). It can be observed that, with an increase in permeability, the non-Darcy factor and the pressure at the entrance of the KCC decreases, while the Forchheimer number and the water inow increases gradually. When the variable n ranges from 10 0 to 10 3 , the Forchheimer number changes from 10 À2 to 10 1 . According to the Forchheimer number and the variable n, the curve can be divided into three phases:   (1) n < 10: the Forchheimer number is smaller than 0.1, which indicates that the inertial resistance is smaller than 10% of the viscous resistance. The effect of the inertial resistance on the ow velocity is weak. The ow in the KCC tends to be Darcy ow.
(2) 10 < n < 580: the Forchheimer number is in the range of 0.1-10, the uid velocity is controlled by both viscous resistance and inertial resistance. However, as the KCC permeability increases, the ow in the KCC is initially dominated by viscous resistance, and then dominated by inertial resistance. When n ¼ 45, the Forchheimer number is about 1, the viscous resistance and the inertial resistance are approximately the same. This phase can be further divided as 10 < n < 45 for weak inertial ow and 45 < n < 580 for strong inertial ow.
(3) n > 580: the Forchheimer number is higher than 10; thus, the inertial resistance in the KCC is at least 10 times that of the viscous resistance. The viscous resistance weakly affects the ow velocity. The uid ow in the KCC is non-linear laminar ow, which tends to be turbulent ow.
The results further illustrate that the Forchheimer equation applied to describe the high non-linear laminar ow can reect the intermediate state of the ow translation from the laminar ow in an aquifer to turbulent ow in a roadway, and also can reveal the essence of the three ow transition in the process of water inrush. When the water pressure in the aquifer is constant, the KCC permeability is an important indicator of transition of the uid ow state in the KCC.

Further discussion
As mentioned above, when the velocity is sufficiently low, the Forchheimer equation can be approximately replaced by the Darcy equation. In other words, the Darcy equation is a particular case of the Forchheimer equation at low velocity. Assuming that the ow in the aquifer can be described by the Forchheimer equation, the critical Forchheimer number, F o ¼ 0.1, is capable of identifying the beginning of non-Darcy ow. 44 Hence, with the increase of the KCC permeability, the ow in the aquifer does not exactly satisfy Darcy's law and the ow in the KCC is also not entirely non-Darcy. Fig. 13 shows the distribution of the boundaries of the non-Darcy effect region under different permeability variation coefficients, n. It can be observed that for ow in the KCC, when n ¼ 1, non-Darcy ow occurred in the top narrow part of the KCC rst, while Darcy ow occurs at the lower wide part. With the KCC permeability increasing, non-Darcy ow expanded to the lower part and the entire KCC is almost non-Darcy ow till n ¼ 100. For ow in the aquifer, when n ¼ 20, there is a small non-Darcy ow region in the aquifer near the KCC. Roughly speaking, the ow in the aquifer satises Darcy's law basically if n < 20. With increasing n, the non-Darcy ow region in the aquifer will extend away from the KCC rapidly. When n increases to 100, the boundary of non-Darcy ow region in the aquifer is about 26 m away from the KCC. With further increase of n, the non-Darcy ow region will grow slowly. This physical process is similar to the spread of drawdown due to a pumping well in hydraulic engineering. 25,26 Overall, during the entire process of water inrush induced by KCC, the ow behavior, both in the aquifer and the KCC, may present Darcy and non-Darcy ow under high water pressure. This primarily depends on the KCC permeability. The higher  the KCC permeability, the more signicant is the non-Darcy behavior and hence, the higher is the risk of water inrush. In engineering practice, the proper ow model should be selected according to the hydraulic conductivity of KCC and the study region of the aquifer to enable a contribution towards a reasonable prediction of water inow and an assessment of potential water inrush. If the study region of the aquifer is wide or the permeability of the KCC is low, the non-Darcy region in the aquifer can be ignored. If the local region of the aquifer is the main study area and the KCC permeability is high, the non-Darcy model can be used to describe the seepage in the aquifer. The discussion can also provide a reference for further establishing a ner model, e.g., the aquifer has a non-Darcy ow boundary. 26

Conclusions
A non-linear ow model was established to simulate the process of water inrush induced by KCC and was used to study the ow behavior of inrush in the Luotuoshan coal mine. The main conclusions can be drawn as follows: (1) To more accurately reect the nature of water inrush induced by KCC in the coal seam oor, a non-linear ow model that couples three elds, i.e., Darcy ow in the aquifer, highvelocity non-linear ow in the KCC and turbulent ow in the roadway, is established. The proposed model reects the essence of transition of ow states during the entire dynamic process of water inrush induced by KCC.
(2) The ow state of the inrush water in the KCC is a crossbreed between Darcy laminar ow, dominated by viscous resistance, and Forchheimer ow, dominated by inertial resistance. The ow velocity increased with time and the Forchheimer ow is initially dominated by viscous resistance, then dominated by both inertial resistance and viscous resistance, and nally dominated by inertial resistance. Therefore, the water inrush induced by KCC is a dynamic ow process with the evolution of ow types.
(3) During water inrush induced by KCC, the lled sediments and broken stone particles in KCC are washed away under the effects of high pressure and high ow rate scouring. Consequently, the permeability of the KCC evidently increases. This further indicates that the Forchheimer equation can reect the ow characteristics over a wide range of permeabilities and ow rates, and can also reect the intermediate state of the ow translation from laminar ow in the aquifer to turbulent ow in the roadway during the water inrush process.

Author contributions
Tianhong Yang, Xiangang Hou and Wenhao Shi contributed to the formulation of overarching research goals and aims; Xiangang Hou and Wenhao Shi analyzed the data; Wenhao Shi and Xiangang Hou wrote the paper; Wenhao Shi and Xiangang Hou contributed to the numerical calculations and simulations.

Conflicts of interest
The authors declare no conicts of interest.