Abstract
The method of two-phase flow is used to analyze the three-dimensional debris flow field via the theory of computational fluid dynamics (CFD). The theory of computational fluid dynamics and finite volume method is introduced briefly, and the rheological characteristic of non-Newtonian fluid is illuminated. The three-dimensional topography model of Niwan gully was established, and the three-dimensional debris flow field was analyzed. The values and distributions of velocity and pressure field of the debris flow 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 flow using the method of two-phase flow. The sediment area of the Niwan gully debris flow can be estimated quickly, and it may have significance for the engineering prevention of debris flow disasters.
1. 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 development of computer technology, the numerical simulation and analysis of non-Newtonian fluids have become possible [8–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–17]. According to the air-liquid 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).
2. 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°2620, E104°4706. 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 km2 (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.
3. 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 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 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., , , , , , , and are six integration nodes, and node is surrounded by the nodes , , , , , and . The letters , , , , , and represent the contact interface of the corresponding volume units. Furthermore, the width of the volume element corresponding to the node in the and directions is and , respectively. The volume integral of equation (2) on the control volume corresponding to 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:
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 is the area of the control volume interface. The generalized source term can be transformed into the following form: where is a constant and is the slope of the function curve at . 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 , , , , , , and are the equation coefficients. The above part is the discretization process of the finite volume method for the governing equations of three-dimensional fluid mechanics problems.
4. 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 (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 and are the progressive viscosity values at extremely low and extremely high shear rates, respectively, is a time-dimensional constant, and is an infinite-dimensional constant. The fluid and the Cross model becomes a power-law model [20]: where is the power-law index, called consistency (unit: ). For the fluid, there is
If we set to zero, we get the form of shear stress
is shear stress between layers when fluid flows, and 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 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).
5. Three-Dimensional Numerical Simulation of Debris Flow
5.1. Modeling and Meshes
The three-dimensional (3D) model of the Niwan debris flow gully is established by the 3D solid modeling method (i.e., from the line to the surface). The modeling codes of AutoCAD and Unigraphics (UG) software build the 3D terrain model (see Figure 4 for model size and mesh division). 3D lines are generated in AutoCAD, imported into UG to generate surfaces, and stretched into 3D entities by stretching commands. The software of ANSYS Workbench was used as the CAD import platform, and the model was divided into the network units in ICEM. The model was 1600 m in width and 3850 m in length and divided into 156,000 tetrahedral grid units.
(a) Mesh
(b) Boundary condition
5.2. Boundary Conditions and Rheological Properties
There are five types of boundary conditions in CFX: inlet, outlet, wall, opening, and symmetry. The ground is set as an adiabatic and rough wall surface, and the value of roughness is set as 0.1, i.e., the average height of the protrusion of the trench bed surface is 0.1 m. The rest of the model is set as an open boundary, and the air inside the calculation area is the average temperature and atmospheric pressure (25°C, Pa). The gravitational acceleration is 9.8 m/s2 (in the negative direction of the vertical axis). Debris flows into the upper inlet of the calculation area. The inlet height is 5 m, and there are two inlet speeds set for debris flow fluid as follows: 1.0 and 3.0 m/s. CFX Expression Language (CEL) in CFX was used to define the viscosity of debris flow fluid. Let the parameter of equation (12) equal 2.0, and the Cross model can be written as follows:
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.
5.3. Calculation Results
CFX mainly reflects the flow field in the form of the streamline. For the 3D terrain model, the velocity distribution contour on the slice plane and isosurface can be displayed. The final calculation results can also be output in tables and curves. The output results after calculation are shown in Figures 5–7. 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 Figures 6 and 7, a part area of the village in front of the Niwan gully will be destroyed by the rush of debris flow if the initial start-up flow speed of the debris flow exceeds 3.0 m/s. Meanwhile, the maximum debris flow velocity in Niwan gully will be 3.0 and 5.8 m/s, if the inlet speed is 1.0 and 3.0 m/s, respectively. The maximum flow velocity of debris flow is located in the primary channel interval of the Niwan gully.
(a) Velocity
(b) Viscosity
(a) Inlet speed equals to 1.0 m/s
(b) Inlet speed equals to 3.0 m/s
(a) Inlet speed equals to 1.0 m/s
(b) Inlet speed equals to 3.0 m/s
6. 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 two-phase 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.
Acknowledgments
The research described in this paper was funded by the National Natural Science Foundation of China (No. 51979225).