Open Access
ARTICLE
Binghui Cui^{1,*}, Liaojun Zhang^{2}
1 College of Civil and Transportation Engineering, Hohai University, Nanjing, 210000, China
2 College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing, 210000, China
* Corresponding Author: Binghui Cui. Email:
(This article belongs to this Special Issue: Advances on Mesh and Dimension Reduction Methods)
Computer Modeling in Engineering & Sciences 2023, 135(1), 357-376. https://doi.org/10.32604/cmes.2022.022309
Received 03 March 2022; Accepted 27 May 2022; Issue published 29 September 2022
The occurrence of landslides is usually introduced by earthquakes or heavy rainfall, and is always accompanied by a large number of casualties and extensive damage [1]. According to the classification of landslides, there is a term “flow-like landslides”, which has the characteristics of fast flow velocity and long run-out distance, that causes severe damage [2–4]. Therefore, great attention should be paid to these flow-like landslides to study the flow mechanism and the sliding characteristics. Due to the vast scale, unpredictability, short duration, and destructive character of landslide disasters, it is quite difficult and dangerous to record the process in real-time. Numerical methods can simulate the flow behavior of such landslides, and can also be applied to the investigation and evaluation of landslides. Most numerical studies on landslides are based on solid mechanics analysis methods, such as the limit equilibrium, finite element method (FEM) [5] and the fracture mechanics method [6]. These methods focus on slope stability analysis, which analyzes the mechanism of landslide damage from the stress-strain perspective. However, these methods are not able to simulate the large deformation state of the landslide, and thus cannot reproduce the whole landslide process.
For physical experiments, it is a very effective and straightforward way of observing landslide characteristics. However, large-scale landslide or debris flow experiments are very expensive and time-consuming. Recently, numerical simulations based on various numerical methods have been widely employed in the simulation of landslides [7]. The discrete element method (DEM) is a feasible method for simulating landslides [8–10]. DEM treats this geotechnical material as a collection of rigid particles governed by Newton's laws of motion. The interactions between particles are represented through appropriate force-displacement relationships that relate to the overlap and contact forces that can occur between particles. However, sufficiently small-time steps are needed to avoid numerical instability, leading to high costs when simulating problems containing a large number of elements. To reduce the computational effort, the grain size in a typical DEM simulation is typically significantly larger than that in a real problem. Another viable approach is discontinuous deformation analysis (DDA) [11,12]. Since DDA was proposed, it has undergone extensive modifications and enhancements to make it more suitable for a wide range of real engineering challenges [13–15]. However, the characterization of fluid flow still needs further improvement [16].
Among the recent advanced numerical methods, SPH is a mesh-free method, a completely Lagrangian method in which the particles carry field variables, (such as mass, pressure, energy, etc.) and move with material velocity [17]. Since it was first introduced in 1977 to solve astrophysical problems in three-dimensional open space [18,19], it has been widely used in different fields due to its meshless and Lagrangian nature. This method provides better robustness and reliability in solving large deformation simulations, owing to the absence of computational grids or meshes [20,21]. Therefore, SPH is more suitable for the simulation of flow-like landslides [22]. The SPH still has some inherent flaws, such as tensile instability, absence of consistency or completeness of interpolation, zero-energy mode, and treatment of boundary conditions, etc. [23]. With the continuous improvement of the SPH method, the method has been successfully applied to geological hazards [2,4,24,25]. Lots of literature [26–28] reported the successful simulation of flow-type landslides as non-Newtonian fluids using the SPH method in combination with the Bingham constitutive model. Most of the existing SPH simulations are based on two-dimensional analysis, which cannot reflect the diffusion and lateral shrinkage, which are different from practical applications. Moreover, the SPH method has recently become more popular, due to the remarkable progress in accuracy, such as the incompressible SPH [29], corrective smoothed particle method (CSPM), and Delta-SPH [30].
Boundary treatment has been the main focus and difficulty of the SPH method, which affects the accuracy and precision of the calculation. Many different approaches have been proposed in the literature to enforce the solid boundary conditions to prevent particle penetration [31–33]. The proposal of a generalized wall boundary condition is an alternative to the dynamic boundary conditions, which can be applied in complex geometries [34,35]. However, the current simulation of non-Newtonian fluids lacks the treatment of fluid-solid boundaries [36]. The literature suggested that the application of free-slip boundary conditions to treat the fluid-solid boundary will make the non-Newtonian fluid unable to obtain resistance from the solid wall boundary in the simulation, thus making the simulation and experiment inconsistent [37].
In this study, a boundary treatment method is proposed that considers friction based on the SPH method. The boundary particle viscous force term is obtained by interpolating the values of the fluid particles near the boundary. In addition, the boundary treatment is verified by comparing the simulation results to the experimental data. The Cross model was employed in the SPH method to simulate the landslides. The rheological parameters were derived from the Bingham model and the Mohr-Coulomb criterion. As a case study, the effectiveness and stability of the application were verified by reproducing the Shuicheng flow-like landslide event in both 2D and 3D that occurred on 23 July 2019, in Guizhou Province, China. The simulation results were compared to the field survey data, showing that the SPH model can provide an accurate analysis of the kinetic characteristics of flow-like landslides, including the flow path, movement velocity, run-out distance and deposition.
For the SPH method, the entire domain is represented as a set of randomly distributed particles with no connection between them. Two main steps are required to obtain the fluid governing equations in SPH form, namely, the kernel approximation and the particle approximation. Firstly, approximation of the field function and its derivative are derived from an integral representation method by using a smoothing kernel function. Secondly, these approximations are replaced by the sum of all neighboring particles in the support domain. It is the kernel approximation and particle discretization that make the SPH method work well for free surface flows and large deformations [17,38].
The integral representation of a function
where
where
There are several kernel functions that are widely utilized. Typically, higher-order kernel functions often enhance the computing accuracy, and substantial time is required. The cubic spline kernel, which belongs to the B-spline family, is one of the kernels that is commonly employed in SPH method [17] and is given by:
where
Due to the density approximation, the pressure field in the SPH formulation for weakly compressible fluids generally exhibits instability and numerical noise. Particularly, in the present case, the noise of the pressure field may be critical. In this study, the delta-SPH method [30] was adopted, which was implemented by adding a diffusive term in the continuity Eq. (5) and it can be written as:
where
where
Regarding the momentum Eq. (6), the total stress tensor
where
In this work, the temperature of landslides is considered as a constant. Therefore, the energy equation is not considered. In the typical SPH method for resolving the compressible flows, the particle motion is governed by the pressure gradient, whereas the particle pressure is determined using an appropriate Equation of State (EoS).
where γ is commonly set to be 7.0;
PySPH is an open-source framework for the simulation of the SPH method [40]. The original PySPH code (version 1.0b1) is able to handle with the free-surface flows, but it is incapable of dealing with non-Newtonian flows. A relevant constitutive model should be developed in the original code to allow for the landslide simulation. In Newtonian fluid dynamics, the shear stress of the fluid can be expressed by the following equation:
where
where
where
Assuming that under the condition of shear rate,
where
Typically,
The Mohr-Coulomb yield criterion with the cohesion
where
Notice that the advantage of the Cross model over the Bingham model is that the effective viscosity is a continuous variable, and the numerical instability is avoided, as shown in Fig. 1.
The boundary treatment has been a main challenge in SPH simulation and affects the calculation accuracy as well as the efficiency. It is vital to select an adequate boundary condition that represents the effect of solid boundary in order to closely describe the dynamic mechanism of landslides, as shown in Fig. 2.
A simple and convenient boundary condition was proposed that can consider the boundary friction. In this work, two types of boundary conditions need to be set. Firstly, when a particle approaches to a solid boundary, the kernel function will be truncated, and an error in the solution result will occur. To address this issue, a generalized wall boundary condition by Adami et al. [46] (see Eq. (23)), was introduced. By extrapolating from the neighboring fluid particles, the pressure on the solid boundary is obtained as:
Secondly, the discrete particles of landslides may slip at the boundary and thus suffer from the slippage resistance (see Eq. (24)). Referring to the concept of pressure interpolation on the solid boundary, the value of the viscosity force on the solid boundary was given by:
Through the above two steps, the no-slip boundary conditions considering friction on the solid boundary was successfully set. When calculating the viscous forces on the boundary, both the dynamic viscosity coefficient can be obtained by interpolation from the fluid, or a fixed viscosity coefficient can be set.
Prediction-correction algorithm was adopted to execute the time integration. For continuity equations and momentum equations of the form, all the values of time-variant quantities are predicted at the time step
Afterwards, these values are updated at another half time step:
At the end of the time step, the values are calculated as follows:
The selection of the magnitude of time step
The flow chart of the SPH with proposed boundary condition is shown in Fig. 3.
3 Verification of the SPH Model
In this section, the proposed boundary condition coupled SPH model is utilized to simulate the granular flow model test conducted by Moriguchi et al. [47] to verify its numerical accuracy. Fig. 4 shows the configuration of the granular flow model. The test was conducted in a tank with the dimensions
Two different flume inclinations were chosen for the two-dimensional simulation, and the parameters are given in Table 1. Each flume inclination is divided into two cases. One case is that the proposed boundary condition is applied to the bottom of the slope. Another case is that the free-slip boundary is applied to all the walls. As the simulation proceeds, a large number of sand particles have a tendency to move from the right side to the left side due to the effect of gravity. Then, the sand particles start to accelerate and hit the measuring wall.
The simulation process with the flume inclination of
To quantitatively analyze the impact of the boundary, the impact force of the flowing sand was measured for different flume inclinations and different boundary conditions. The simulated time histories of the impact force of the flowing sand are presented in Fig. 6. Comparing the results with the tested data, these two figures show that the peak impact force value of impact forces can be well captured by the SPH model. As excepted, comparing different boundary cases, the peak impact force generally agrees well with the tested data when the friction boundary is considered.
4.1 Geological Setting and Flow-Like Landslides Features
The Shuicheng flow-like landslide event was triggered by intense rainfall in July 2019 (Fig. 7a). A total of three heavy rainfall events occurred in the week before the occurrence of landslides, with a cumulative rainfall of 288.9 mm. Under the influence of continuous heavy rainfall, the landslide began to lose stability and rapidly slid along the slope at high speed. After roughly 700 m of landslide movement, two diversions were formed after hitting a small ridge and subsequently impacting and burying the residential buildings on the lower part of the slope. The mass flow of landslide spread along the slopes for about 1490 m and was deposited in the gully depressions. Due to the suddenness of the landslide and the lack of protective measures such as blocking, the disaster destroyed 21 houses and 43 people were killed.
The Guizhou Provincial Geological Disaster Emergency Technical Guidance Center immediately carried out the field investigations following the occurrence of landslide. The primary characteristics of the flow-like landslides are depicted in Fig. 7b as well as a typical cross-section of the Shuicheng gully in Fig. 5c. According to the field investigations, the entire landslide is approximately 1300 m long, with a large distance of around 460 m between the front and rear edges, and a volume of about 191.2 × 104 m3. According to Zhou et al. [48], the average density of the Shuicheng landslide mass was about 2000 kg/m3. Depending on the test results, the angle of internal friction and the cohesion of the flow-like landslide can be set approximately as 30° and 30 kPa, respectively. In the previous simulations, the Bingham model was widely applied to simulate the landslides or debris flows considering a range of dynamic viscosities from 20 to 500
where
According to the simulation results, the whole process duration was about 90 s. Fig. 8 presents the slope configurations at different time points, reflecting the movement of landslides. The particles slipped down from the top of the landslide and then moved along the steep slope by gravity. Finally, these particles reached an equilibrium state and were deposited in the lowlands. The overall performance of the sliding was accelerated motion in the period 0–30 s and decelerated motion after 30 s.
Landslides usually have a long distance, and the measurement and selection of various parameters may vary greatly in different locations. Therefore, it is meaningful to conduct a sensitivity analysis to determine how various rheological parameters affect the simulation results. Table 2 lists the eight calculation cases with different rheological parameters. To explicitly quantify the differences of the sensitivity of different parameters
where
Since the 2D simulation cannot reflect the diffusion and lateral contraction, 3D simulation on the real complex terrain is necessary. In this section, the 3D terrain is generated from the original topographic map at a scale of 1/2000 and homogenizes the mesh. The previously created mesh is then transformed into a sequence of boundary particle particles. Similarly, the landslide is discretized into a succession of particles each with its own set of properties. In this simulation, the particle distance is set to 7.5 m, resulting in 12825 boundary particles and 7010 fluid particles. The strength parameters adopted in the 3D simulation are the same as in the 2D simulation (Table 3). Based on this model, Shuicheng landslide was numerically simulated in 3D terrain, and the results are shown in Fig. 10. The color of the particles in the figures represents the sliding velocity at various instants.
The results suggest that the front flow takes about 28 s to travel nearly 520 m to reach the isolated island. The front flow velocity reaches a maximum with an average velocity of 29.8 m/s. Due to the good stability of the soil near the isolated island, it diverts the debris to the sides and forms a divergence. The diversion of the stable slope protects the lower residential houses and creates an “island of safety”. The existence of isolated islands has dramatically slowed down the landslides. Therefore, the front flow remains stable at about 90 s and is deposited in the lowlands. Eventually, the estimated total volume of deposition is 1.1 × 106 m3, which is less than the field investigation value 1.91 × 106 m3. The reason might be that the simulation process does not take into account the entrapment process. In addition, the simulated landslide paths and deposited areas matched well with the results of field investigations. The topographic data, in this case, were obtained from the 1/2000 scale topographic map, and it is considered that more accurate topographic data would help to obtain better simulation results.
The displacement and velocity time histories of the landslide front and rear were investigated to gain a better understanding of the landslide characteristics in the present simulation. Fig. 11 depicts the velocity time series of the landslide's front and rear. At around 20 s and around 45.5 s after failure, the flow-like landslide front and rear reach a maximum velocity of 29.2 and 26.5 m/s, respectively. Afterwards, an isolated island was encountered, the landslide diverged and the path was changed, and the speed gradually decreased. After moving for 70 s, the landslide came to a halt and was deposited in a lower area. It should be noted that the changes in velocity at various moments represent the complexity of the 3D terrain. Fig. 12 shows the maximum displacements of the front and rear during the motion were 1099 and 859 m, respectively, at the time of 90 s after failure. Ma et al. [51] conducted extensive field investigations and calculated the maximum velocity of 25 m/s2, and a landslide duration of 60 s using the empirical equations. Both values turned out to be slightly smaller than the simulated results. Quantitative validation of the landslide simulation based on the field measured data and empirical equations were challenging. The error may be triggered by a combination of natural and artificial reasons, which makes the simulation more challenging. In addition, it can be seen that the flow's pressure field is formed smoothly. Fig. 13 depicts the evolution of the flow-like landslide pressure field over time. As previously stated, delta-SPH successfully corrects pressure field noise, particularly for the high viscosity fluids.
The landslide material is discretized into a succession of SPH particles of almost the same diameters in this model. The greater the number of particles and the narrower the particle spacing, the greater the computing time required. It should be mentioned that, in the current code, the SPH module is accelerated based on Open Multi-Processing (OpenMP). In theory, the higher the number of particles in the SPH, the more accurate the simulation will be, but a balance needs to be reached with the needed calculation time. Fig. 14 depicts the relationship between the program runtime and thread count. The results show that the computational efficiency of the proposed SPH model improves as the number of threads increases. Although the OpenMP code is more computationally efficient than the serial code on the CPU, the computation time consumption is large when more particles exist.
A rational simulation of flow-like landslide propagation process contributes to hazard analysis and disaster mitigation design. In this study, the Navier-Stokes equation combined with the non-Newtonian fluid rheology model was utilized to investigate the dynamic behavior of the flow-like landslide.
1. A simple and convenient boundary condition was proposed that can consider boundary friction. By using the proposed boundary condition, the simulation results are in good agreement with the experimental results. Moreover, since this operation is limited to boundary particles, it will not increase the computational overhead.
2. The Shuicheng flow-like landslide, which occurred in Guizhou in 2019, was studied as a case. The parameter sensitivity is analyzed under the 2D model, and the results show the shearing strength parameters have more influence on the computing accuracy than the coefficient of viscosity.
3. A 3D model of Shuicheng landslide was constructed and computed that corresponded to the site conditions. In terms of the solid volume and deposition area, the computed results are generally consistent with the field investigations. In addition, the characteristics of the front and the rear flow velocity and zone of influence of the flow-like landslide were analyzed, which will help the mitigation or countermeasure design work and offer evidence for the hazard assessment.
4. The approach provides an effective tool for studying the dynamic behaviors of landslides. However, the current study does not take into account the entrainment process of flow-like landslides. Large-scale practical problems usually require massive particles, and future research should be directed to GPU-accelerated modules.
Acknowledgement: The authors appreciate the valuable discussions with Zhang Weijie (College of Civil and Transportation Engineering, Hohai University, China) and Shi Wenbin (College of Civil Enginnering, Guizhou University, China) on topics of landslides.
Funding Statement: The authors received no specific funding for this study. after this line.
Conflicts of Interest: We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.
References
View
Download
Like
Modeling and Simulation of Sealing Spray Application Using Smoothed Particle Hydrodynamics
Robert Rundqvist, Andreas Mark,...Non-Newtonian Blood Flow in Left Coronary Arteries with Varying Stenosis: A Comparative Study
PoojaJhunjhunwala, P.M. Padole,...CFD Analysis of Pulsatile Flow and Non-Newtonian Behavior of Blood in Arteries
P. Jhunjhunwala, P.M. Padole,...Cytoplasmic Motion Induced by Cytoskeleton Stretching and Its Effect on Cell Mechanics
T. Zhang