Dynamic Simulation and Analysis of Large-Scale Debris Flow Field

The method of two-phase ﬂ ow is used to analyze the three-dimensional debris ﬂ ow ﬁ eld via the theory of computational ﬂ uid dynamics (CFD). The theory of computational ﬂ uid dynamics and ﬁ nite volume method is introduced brie ﬂ y, and the rheological characteristic of non-Newtonian ﬂ uid is illuminated. The three-dimensional topography model of Niwan gully was established, and the three-dimensional debris ﬂ ow ﬁ eld was analyzed. The values and distributions of velocity and pressure ﬁ eld of the debris ﬂ ow have been obtained; the relationship between the shear strain rate, viscosity, and velocity was analyzed. The results suggest that it is more accurate to catch the free surface of debris ﬂ ow using the method of two-phase ﬂ ow. The sediment area of the Niwan gully debris ﬂ ow can be estimated quickly, and it may have signi ﬁ cance for the engineering prevention of debris ﬂ ow disasters.


Introduction
Geological disasters such as debris flow and landslides frequently occurred recently due to global climate change [1]. As a geological disaster, debris flow has been included in the "International Decade of Natural Disaster Mitigation" project by the United Nations due to its sudden outbreak and incredible destructive power [2]. China is one of the countries with maximum events of debris flow disasters. The worst debris flow disaster in China happened in Zhouqu, Gansu province, on Aug. 8,2008, and it killed thousands of people and affected dozens of administrative villages [3]. Debris flow fluid belongs to non-Newtonian fluid with significant rheological characteristics [4]. The granular flow expandable model was proposed by Takahashi to describe debris flow movement [5]. Chun proposed a non-Newtonian plastic expansive body model that summarizes the dynamic models of debris flow [6]. Cui et al. proposed a more practical fluid movement model of debris flow [7]. Due to the progress of fluid mechanics, computational fluid dynamics, and rheology and the rapid develop-ment of computer technology, the numerical simulation and analysis of non-Newtonian fluids have become possible [8][9][10]. Choi et al. analyzed the three-dimensional debris flow field and its influence on the gully surface via the smoothed particle hydrodynamics (SPH) method [11]. Literature [12] calculated and analyzed the flow and deposition of debris flow based on FLO-2D simulation. It is of great significance to process the numerical simulation and analysis of debris flow via the theory of computational fluid dynamics (CFD) [13,14]. The variation relationship between velocity, pressure field, and terrain of the debris flow gully is significant to study further the dynamic movement mechanism and two-phase characteristic of debris flow fluid [15][16][17]. According to the airliquid two-phase flow method, the free surface of fluid can be captured more accurately [15]. The air-liquid two-phase flow method is used to track the free surface of debris flow fluid to determine the silting range of debris flow more efficiently. In this paper, a case of the debris flow gully is used to analyze the flow field of debris flow via the theory of computational fluid dynamics and finite volume method (FVM).

Background
The debris flow gully on the left bank of Bailong River is selected as a typical place for simulation analysis. Niwan gully is located in Houba Village of Longnan City. The coordinates of the gully are as follows: N33°26 ′ 20 ″ , E104°47 ′ 06 ″ . The drainage basin has a maximum elevation of 2516 m and flows into Bailong River at an altitude of 1037 m with a relative elevation of 1479 m. The basin has a steep topography on both sides of the gully and covers about 11.53 km 2 (shown in Figure 1). The topography of the Niwan gully is a typical medium and low mountain canyon landform of tectonic erosion and denudation. The main gully is 7.66 km long and has an average longitudinal dip of 193.2‰. The surface of the gully covers thin loess, and the soil erosion in the gully is severe. Niwan gully has seasonal water flow, and the heavy rain (mainly concentrated from May to Sep.) is the primary source of hydrodynamics. According to the data, the minimum rainfall intensity of debris flow in the mountain area of Longnan is 15-20 mm/h, and the maximum precipitation in this area is 40 mm per hour. Rainfall erodes the slope surface and bank of the gully, and the debris flow may easily be triggered.

Theory and Method
This paper uses computational fluid dynamics (CFD) software CFX to calculate and analyze the characteristics of the debris flow field. CFD's basic principle theory uses a series of discrete numerical methods to approximate the continuous physical field via the governing equations of fluid flow (conservation equations of mass, momentum, and energy) [14]. The finite volume method is the main discretization method in computational fluid dynamics (CFD), unlike the finite element method. The finite volume method divides the analysis object into several control volumes, integrates the governing equation of the problem to each control volume, and approximates the physical quantity of the volume    2 Geofluids interface by the physical quantity on the node through the interpolation method. The three basic governing equations in fluid mechanics can be written as follows [14]: where ф is the unknown variables (speed and temperature), Γ is diffusion coefficient, and S is the generalized source term. Equation (1) can be written as components: Figure 2 is a schematic diagram of volume unit division in the computational domain. Each volume element contains an integration node, i.e., P, W, N, E, M, B, and T are six integration nodes, and node P is surrounded by the nodes W, N, E, M, B, and T. The letters w, n, e, m, b, and t represent the contact interface of the corresponding volume units. Furthermore, the width of the volume element corresponding to the node P in the x and y directions is Δ x and Δy, respectively. The volume integral of equation (2) on the control volume corresponding to P is The finite volume method mainly adopts the central difference method to discretize the parameters. That is, the linear interpolation of the variables in equation (2) on each control volume interface is as follows: Similar treatment should be made for other contact surfaces. Substitute the above interpolation form into equation (3), then the transient term can be written as follows: The Gauss divergence theorem is used to transform the volume fraction into a surface integral, where the convective term can be written as follows: Newtonian behavior Non-Newtonian (shear thinning) behavior

Geofluids
Substitute the interpolation into equation (6): The diffusion term related to the gradient can be obtained as follows: Substitute the interpolation into equation (6) which can be obtained as follows: where A is the area of the control volume interface. The generalized source term can be transformed into the following form: where S C is a constant and S P is the slope of the function S curve at P. Then, the generalized source term is as follows: After substituting the discrete expressions of the above items into equation (3) and adding the constant terms, a set of algebraic equations can be obtained as follows: where a P , a W , a E , a S , a N , a B , and a T are the equation coefficients. The above part is the discretization process of the finite volume method for the governing equations of threedimensional fluid mechanics problems.

Rheological Properties of Non-Newtonian Fluids
The rheological properties of non-Newtonian fluids refer to the relationship between shear stress and shear rate when they are subjected to shear deformation after flow [18,19].
In the same way that viscosity is defined for Newtonian fluids, viscosity for non-Newtonian fluids can be expressed as a function of the shear rate du/dy (also known as the velocity gradient): In recent years, many non-Newtonian fluid viscosity formulas have been proposed by scholars of various countries, among which the Cross model has strong applicability [20]: where μ 0 and μ ∞ are the progressive viscosity values at extremely low and extremely high shear rates, respectively, K is a time-dimensional constant, and m is an infinitedimensional constant. The fluid and the Cross model becomes a power-law model [20]: where n is the power-law index, called consistency (unit: Pa × s n ). For the fluid, there is If we set n to zero, we get the form of shear stress τ is shear stress between layers when fluid flows, and τ y is yield shear stress. This formula is the Bingham model. The viscosity of the non-Newtonian fluid with shear rate can be divided into shear thinning and shear thickening. Shear thinning is a flow in which viscosity decreases with the increase of shear stress or shear rate and vice versa. Caricchi et al. proposed the magmatic rheological tests of different crystal components showing that the fluids turn to Newtonian fluid characteristics when the shear rates are lower than 10 -5 . When the shear rate is between 10 -5 and 10 -3 , it shows non-Newtonian fluid 3 6 Geofluids characteristics, and the viscosity decreases with the increase of the shear rate. The shear-thinning phenomenon of pseudoplastic fluid occurs. The fluid shows Bingham fluid characteristics when the shear rate value is great than 10 -3 [21], and the curve of magmatic rheological characteristics is shown in Figure 3. At high shear rates, viscosity approaches a constant, in which case the power-law model is no longer applicable. The Cross model can reflect the rheological characteristics of non-Newtonian fluids in a shear rate of 4 or 5 orders of magnitude at a high shear rate [21]. The Cross model is used to describe the rheological characteristics of debris flow fluid. The viscosity of debris flow fluid follows equation (12).

Three-Dimensional
Equation (18) is the viscosity expression of debris flow fluid in numerical simulation. Relevant air parameters contained in the software are directly used in the calculation, while debris flow fluid is added as a customized fluid. The calculation parameters and values of the Cross model are shown in Table 1.  Figure 5 is the distribution contour of velocity on the free surface between debris flow fluid and air calculated by the air-liquid layered two-phase flow method. Figure 6 shows the velocity contour on the free surface with the satellite image of Niwan gully. Figure 7 shows the velocity contour on the free surface with a remote sensing image of Niwan gully. The inlet speed for the analyses can be regarded as the debris flow's initial start-up flow speed. The initial speed of the debris flow is influenced by the surface runoff caused by rainfall. The free surface contour of the debris flow field combined with the satellite image of the debris flow gully and the submerged area of the debris flow disaster will be determined and evaluated. From

Conclusions
A three-dimensional debris flow field with a free surface is analyzed via the computational fluid dynamics and finite volume methods (FVM). The air-liquid two-phase flow method is used to capture the free surface of debris flow and determine the siltation range of the three-dimensional debris flow field. The conclusions are as follows: (1) The free surface of the three-dimensional debris flow field can be captured accurately by the layered twophase flow method. The distribution of velocity and viscosity on the free surface can be obtained efficiently. Therefore, it can provide a basis for judging debris flow's siltation range and preventing and controlling disasters (2) The relationship between the flow field of debris flow and the change of gully topography can be obtained through calculation and analysis of FVM, which has theoretical significance for the depth study of the rheological characteristics of non-Newtonian fluid (3) Various non-Newtonian fluid rheological models can be customized with the computational fluid dynamics theory. The influence of the roughness of the gully surface can be considered, and more research methods are proposed to study debris flow movement and dynamics (4) The influence of solid material with large diameters (such as gravel) on debris flow cannot be predicted by FVM with the theory of computational fluid dynamics. This work can be further analyzed using the discrete element method

Data Availability
Some or all data, models, or code generated or used during the study are available from the corresponding author by request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.