This study proposes a solid-fluid-coupled Smoothed Particle Hydrodynamics model to investigate the behavior of pipe in the liquefied sand after failure. In this model, the liquefied soil is simulated as a Bingham fluid material combined with the equivalent Newtonian viscosity. Then the pipe culvert is treated as an elastic solid material, and the culvert-soil interaction is assumed as the solid-fluid coupling force. Verification has been conducted through a dam-break model to demonstrate the accuracy of simulating the flow behavior of liquefied sand. At the same time, the applicability of the SPH simulation in the two-phase-coupled force has been checked compared with previous researches. After that, parametric studies are performed on the lateral force and movement of pipe culverts in the flowing sand. The derived results demonstrate that the culvert size and buried depth are the significant factors when estimating the maximum lateral force of culvert-soil system, while the slope angle severely affects the run-out distance of culvert and the velocity of flowing sand.

1. Introduction

Underground pipelines or culverts suffered great destructions when subjected to the lateral deformation of liquefied soil [14]. In practice, culverts are composed of attached segments of pipe or cast-in-place monolithic structures but are shorter in length than pipelines. There have been several cases of severe damage primarily caused by the liquefaction-induced embankment penetration or spreading. Among them, the overflow conduit failure at Jensen Water Filtration Plant during the 1971 San Fernando earthquake [1] and the culvert failure near Lower San Fernando Dam during the 1994 Northridge earthquake [2] represented typical phenomena of lateral deformation-linked damage. The first case was that the concrete pipe culvert with a diameter of 2.44 m was pulled apart at the joints and suffered lateral displacements as great as 1.8 m measured at the face of compacted fill near the conduit outlet. And the second case was that the corrugated metal pipe culvert with a same diameter was separated at pipe joints and segments were driven with a large lateral distance. Clearly, the culverts have insufficient strength in joints to fully resist these ground displacements and would lose lateral restraint when suffering large deformation of soil once these joints are damaged. However, it is difficult to use experiments or field investigations to detect these destructions and behaviors after failure in detail due to their concealment feature. Therefore, it is much necessary to conduct a numerical analysis on this kind of phenomenon and provide disaster prediction, prevention, and mitigation for the underground engineering.

Numerous experimental or numerical studies have been performed to evaluate the culvert-soil systems [58]. Katona [6] proposed a plane-strain finite-element method to evaluate the failure performance of buried culverts under seismic loading. The dynamic centrifuge model tests had also been conducted to investigate the seismic behavior of culverts, while FLAC program was used to investigate the soil properties and configuration effect on the mechanical behavior of culverts [8]. These studies mainly focus on the static or dynamic behavior of culverts under normal operation and have not been involving the postfailure behavior of culverts in the flowing soil.

From previous researches, it could be found that the behavior of liquefied sand is very close to the fluid as the pore water pressure occupies the effective stress on the soil skeleton [9, 10]. Thus it has been validated that the Bingham model could represent the large-deformation behavior of liquefied sand in the numerical analysis [11, 12]. The grid-based methods like FEM and FDM would encounter some numerical difficulties (e.g., severe mesh winding, twisting and distortion), when it comes to large-deformation problems. As a mesh-free particle method based on the Lagrangian description, the SPH method has an advantage in the simulation of large deformation [13] and has been widely used in the coupling modelling of solid phase and fluid phase [1418]. In the aspect of solid-fluid coupling, a soil-water-coupled SPH modeling method had been developed to analyze the flowing behavior of liquefied soils, where liquefied sand was modeled as the water-saturated porous media [14]. At the same time, this model had been used to simulate the debris flow and regarded fluid-structure interaction as the repulsive solid-fluid force to assess the impact force of debris flows [15]. Meanwhile a two-phase-coupled SPH model in coordination with a constitutive model for unsaturated soils had been proposed to conduct the unsaturated seepage analysis of dike [16]. On the whole, the SPH method has been applied widely to the simulation of liquefied soils in existing literatures. Yet, little study has been undertaken on the interaction between structures and flowing soil since it is difficult to compute the flowing soil pressure on the structure.

The aim of this study is to propose a SPH model for the culvert-soil system and to calculate the lateral force and movement of culvert without restraint from the third direction (i.e., unconfined culvert) in liquefied sand. In the SPH model, the culvert is treated as an elastic material while the flowing sand as the Bingham fluid, where equivalent Newtonian viscosity is used. The fluid pressure of liquefied soil on culverts is added into the momentum equation for considering the interface interaction. Then, a computational algorithm is used to solve the contact penetrating problem. The flow characteristics of liquefied soil are the key point in the whole simulation process; thus the accuracy of the proposed SPH model on the fluid phase is verified by simulating the dam-break problem. At the same time, the solid-fluid coupling is validated by comparing the impact force of a flowing test with measured data. Then, this model is applied to the lateral force and movement investigation of the unconfined culvert-soil system, to reveal the effect of culvert size, buried depth, model slope, and boundary condition.

2. SPH Numerical Framework

2.1. Basic Concepts

In the SPH method, the problem domain is discretized by a set of unconnected particles. These particles carry field variables such as density, stress, and velocity. Through SPH, the field function in continuous form is represented using an interpolation function and the function value of each particle is approximated using averages of function values at all neighboring particles as follows:where r is the position vector, N is the number of supporting particles which influenced the computing particle, f(r) is the function value, m is its mass, ρ is its density, Wij = W (ri - rj, h) is a smoothing kernel function, and h is the smoothing length which determined the number of supporting particles. The subscripts i and j present the ith and jth particle in this study.

Among the kernel functions of previous studies [1921], the cubic spline kernel function proposed by Yang and Liu [22], which has reasonable computational cost and could eliminate the effect of stress instability when computational domain was in compression, is used in this simulation as follows:where s = /h, and h is the smoothing length. The smoothing length controls the size of search area and number of supporting particles around particle i. In this paper, the smoothing length is twice the particle spacing.

In addition, the link-list searching algorithm to search for the neighboring supporting particles is used in the proposed SPH model.

2.2. Governing Equations

Figure 1 gives the schematic diagram of the proposed SPH model. The culvert is computed as a solid material. The governing equations for traditional SPH include the continuity equation and the momentum equation as follows:where is the stress tensor, v is the velocity vector, b is the vector acceleration for gravity, and ▽ denotes the Nabla operator.

These could be approximated by applying (4) and (5) to (1) and (2),

To improve the overall stability of the simulations, the artificial viscosity is introduced in the above momentum equation,where is Kronecker’s delta and Пij is an artificial viscosity used to mitigate the numerical oscillation and avoid the penetration between particles [23]. Пij could be obtained from the following equation,where α is 0.001, β is 0, rij = ri - rj is position vector increment, and cs is the sound speed.

In the aforementioned classical SPH theory, the application for solid mechanics has some problems like boundary deficiency, particle-distribution sensitivity, and ill-conditioned matrices. These problems are caused by boundary truncation of support domain, discrepancy between kernel and particle approximation, as well as highly disordered particle distribution. The insufficient particles around the interface between the culvert and the sand lead to the low accuracy in the SPH simulation. To address these problems, the normalized particle approximation referring to the decoupled finite particle method (DFPM) proposed by Zhang and Liu [24] is adopted in this paper. The corrected form of particle approximation could be obtained by the following formula:

To moderate the velocity of the computing particle to match the average velocity of surrounding particles, the ‘XSPH’ approach [25] is used in this context to revise the particle velocity as follows:where vi is the revised velocity vector of particle i, vi’ is its original velocity vector, and ε is a constant between 0 and 1. The value of ε is a representation of how close the revised velocity to the average velocity. α is 0.5 for culvert and 0.3 for soil in this work. And the choice of these values could avoid one particle penetrating into adjacent particles in high-speed flows.

2.3. Constitutive Model and Contact Algorithm

In this work, the constitutive model of culvert is an elastic model as follows:where is the stress tensor, is the strain tensor, I1 is the first strain invariant, λ is the Lame constant, and G is the shear modulus.

The particle approximation of strain rate tensor for the culvert could be written as

And the stress rate tensor is calculated as follows:where is the Jaumann stress rate tensor calculated from the constitutive model and is the spin tensor, obtained from the following:

The flowing sand is regarded as a kind of quasi-incompressible fluid with shear strength in this simulation. The following Bingham model, which is used for the landslides simulation [26], is adopted to describe the behavior of flowing sand in this paper,where is the stress tensor, is the strain rate tensor calculated by (19), p0 is the hydrostatic pressure, and is the equivalent Newtonian viscosity [27] shown by following:in which τr is the residual strength, cr is the residual cohesion, φr is the residual internal friction angle, η0 is the viscosity, is the second invariant of the deviatoric strain rate, and p is calculated from the following state equation [19]:where γ is 7.0. The lower this value γ is, the higher the compressibility of water is. In the above equation, even small changes in density bring about large changes in pressure. Actual calculations using this value (γ = 7.0) illustrate that changes of density were less than 0.1% [28], which means a good approximation of incompressibility.

In the aforementioned part, the liquefied soil could be treated as a kind of material like fluid. To some extent, the interaction between culverts and liquefied soil could be investigated as coupling between the deformation of culverts and the motion of liquefied soil. In the past, the interface interaction is modeled as a surface contact in pipe-soil systems [2931]. The pressure of liquefied soil, as a kind of fluid pressure, on structures is seldom considered in most liquefied soil-structure interaction analyses.

In saturated soil, the total stress acting on soil skeleton is calculated from two parameters (u, pore water pressure and σ’, effective stress) [32]. The same is true for structures in liquefied sand. Concerning the momentum equation involved by Zienkiewicz and Shiomi in 1984 [33], (8) could be presented as follows:where the stress is expressed as follows:

The interaction term Ri could be expressed aswhere n is porosity, k is permeability coefficient, and the gravitational constant g is set equal to 9.8 m/s2. The value of Ri is set to 0 when the sand particles are beyond the supporting area of the culvert particles. The contact penetration is also one of the key problems in this simulation. To prevent the soil particles penetrating through the pores among the culvert particles, the soil particles around the culvert should be at a somewhat smaller acceleration. Actually, when the soil particles around the interface area tend to approach the culvert particles, the acceleration of the soil particles could be revised. The following formula is used to judge the tendency:

In this paper, the revised acceleration of soil particles is adopted as follows while the value of Formula (25) is negative:where ai’ is the particle acceleration updated from (22) and β is a constant parameter related to the decreasing range of acceleration. This is an improvement to the particle penetrating problem. The soil particles could not embed into the culvert by using the interaction force algorithm and the acceleration correction method in this simulation.

2.4. Boundary Condition

When boundary particles are in the influencing domain of nonboundary particles, the nonboundary particle velocity moving to the boundary particle should be decreased in order to avoid the penetrating of nonboundary particles. The boundary treatment method used in this work is similar to that used in modeling low Reynolds number incompressible flows [34]. This boundary treatment makes the numerical simulation result closer to the actual response of the liquefied soil, so that the velocities of soil particles around boundary particles became gradually slower. The boundary effect should not be considered if the following inequality is satisfied:in which vi is the velocity vector of nonboundary particle i relative to neighboring boundary particle, and xi is the equivalent displacement vector.

If relation (27) is not satisfied, the velocity of nonboundary particle i should be corrected as follows:where vi is the corrected velocity vector, is the previous velocity vector, and β is calculated from the following equation,in which is 1.6, a is 0.3, b is 0.1, and d is the initial particle space. As shown in Figure 1, the corrected velocity of particle B with one layer of particles, , has higher vertical absolute component than the corrected velocity with two layers, . In other words, the latter is more effective than the former in responding to the velocity correction. Thus, the two layers of boundary particles are built in this paper.

3. Verification of the Proposed Model

As mentioned before, the liquefied sand could be simulated as a Bingham fluid material, and the culvert-soil system could be treated as the solid-fluid system. Definitely, the fluid phase and the coupling force are the key points to verify the reliability of the proposed SPH model. Here the authors present the verification of fluid phase and coupling force using dam-break tests and flowing sand tests.

A dam-breaking case in a smooth rectangular tank was simulated by using the proposed SPH model. The diagram of the dam-break case is illustrated in Figure 2. For simplicity, the nomenclatures are the same as parameters given in above equations. The model tank was 307 cm long and 171 cm high. Two SPH simulations were conducted: one was 5.4 cm in width and 5.4 cm in height and another was 2.9 cm in width and 5.8 cm in height. The fluid density was set at 1×103 kg/m3, residual strength indexes were cr = 0, φr = 5°, and the equivalent Newtonian viscosity η was 1.0×10−3 Pa·s. The fixed computational step was Δt = 1×10−4 s with an initial particle spacing of 0.1 cm. At the beginning stage, the total numbers of particles in the simulation were 5296 and 3701 while the numbers of moving particles were 3364 and 1770, since the rest of the particles corresponded to the boundaries of the tank. Simulated results are compared with the reported testing data [35], as shown in Figure 3. It could be seen that the distances of the free surface using the proposed SPH model are predicted well, thereby proving that the SPH model has relatively high reliability to the flowing simulation of fluid phase, which is important to the large-deformation problem of culvert-soil system.

To validate the interaction between fluid particles and solid particles, the flowing soil test was simulated to compare with the test results [27]. The flowing soil was simulated as Bingham fluid. The particle number of soil particles was 1612, the density was 1.8×103 kg/m3, and the equivalent Newtonian viscosity η was 1.0 Pa·s. Meanwhile, the calculation time step was 8.0×10−5 s and total steps were 2.5×104. The sketch of flowing soil model was as Figure 4, in which a pressure gauge consisted of solid particles was placed in the end of the SPH model box. After that, time series of impact force were measured for slope angles with 45°, 50°, 55°, and 60°. In the simulation, fluid particles were at rest originally and could flow under the action of gravity once the calculation began. Figure 5 gives the comparison of tested and simulated stable and maximum impact force on solid particles. It could be seen that both the stable and the peak impact force of SPH simulations had a good agreement with other research. This indicates that the proposed SPH model has good accuracy in calculating the interaction force between fluid material and solid material.

4. Model of the Culvert-Soil Systems

4.1. Introduction of the Simulating System

In this part, the culvert was defined as a hollow and segmented tubular structure separated at joints while the surrounding soil was regarded as liquefied soil. Thus, this investigation did not consider the part that the force tries to grab and hold the culvert in place. Although this assumption did not satisfy the widely accepted practice, simulated conditions had really happened [1, 2]. Thus this study could provide meaningful conclusions for the engineering practice to some extent.

The geometries of SPH cases were shown in Figure 6, and the flowing depositions were indicated with dash-dot line. In Figure 6, the coordinate system was x axial in the horizontal direction and y axial in the vertical direction, with an incline angle θ. For the case in Figure 6(a), the right wall was far away from the sand profile, while the right wall was close to the sand profile for case in Figure 6(b).

Parameter values in the simulation are listed in Table 1, where the density, diameter, and thickness of pipe culvert are referring to the guidelines of International Tunnel Association [36], in order to ensure that the results could be used for the practice of large pipelines or tunnels. Meanwhile, the particle density of culverts was calculated by using the principle of equivalence of mass [37]:where ρc = 2 650 kg/m3 was the density of culvert referring to the Guidelines [36], tc was the equivalent thickness, ρi was the density of particle i, and di was the particle space. Here, tc was 0.25 m, ρi was 1840 kg/m3 for D = 3 m, while tc was 0.45 m and ρi was 2626 kg/m3 for D = 6 m. And the particle density of sand was 1640 kg/m3.

In this study, the width and height of columns were denoted as L0 and H0, while the diameter and buried depth of pipe culvert were denoted as D and H, respectively. For the case in Figure 6(a), the width of model box quintupled the height H0 of sand profile and the computation domain size was L0 = 30 m and H0 = 60 m. As shown in Table 2, two categories were simulated: one was 3 m in diameter with 20027 sand particles and 44 culvert particles; another had a diameter of 6 m with 19791 soil particles and 112 culvert particles. The computational time step was Δt = 1.25×10−4 s and initial particle space was 0.3 m for both categories in Table 2. In addition, the case in Figure 6(b) had the same simulating categories and schedules as the case of Figure 6(a), also presented in Table 2.

4.2. Results and Discussion

Figure 7 illustrates behaviors of culvert-soil systems under different diameters, which include the changes of pipe culvert center position, lateral forces on culvert, and lateral and vertical average velocities of flowing sand over time. As shown in Figure 7(a), for the case of H/D = 0.5 and θ = 0°, diameter is the comparable variable and the flowing frontier is far away from the right boundary (regarded as without right boundary in this simulation). It is shown that the run-out process of culvert with a large diameter (D = 6 m) motivated by the sand flowing is similar to that of a small diameter (D = 3 m) and their culvert centers are nearly overlapped in the last (t = 10 s). Besides, Figure 7(b) shows that the magnitude of lateral force on the culvert of a larger diameter is considerably larger (nearly ten times) than that of a smaller diameter at the initial stage (t < 3 s) after flowing. Figures 7(c) and 7(d) compare the diameters effect on the sand flowing velocities along with velocities of the simulation without culvert. It could be seen that the pipe culvert diameter has almost no effect on the sand flow velocities and the time of maximum velocities in case without the culvert is slightly later than that with culvert.

Figure 8 presents the effect of the buried depth on culvert-soil interactions. From Figure 8(a), moving curves of the culvert center with different buried depths are almost parallel till that at the time (t > 7.5 s) the culvert in a large buried depth (H = 7.5 m) gradually approaches the culvert in a small buried depth (H = 1.5 m, 4.5 m). Figure 8(b) shows that lateral forces on the culvert increase with rising buried depths at the initial stage (t < 3 s). It could be seen from Figures 8(c) and 8(d) that soil flowing velocities are basically unchanged for various buried depths and the time of maximum velocities in case without culvert is slightly later than that with culvert.

The effects of slope angles on the culvert-soil systems are shown in Figure 9. From Figure 9(a), larger slope angle brings about further culvert flows within the same time. Figure 9(b) indicates that the lateral force–time curves under different slope angles are close. The major difference between them is that the appearance of maximum lateral force over time is distinct. For example, the maximum lateral force on the culvert of a large slope angle (θ = 20°) is at the beginning of flowing. However, the lateral force of small slope angles (θ = 0°, 5°) increased firstly and then decreased, while the maximum lateral force is at t = 1.25 s. It is observed from Figures 9(c) and 9(d) that the lateral velocity and vertical velocity have a direct relationship with the slope angle, which is that increasing slope angle means increasing velocities, especially for the lateral velocities. Meanwhile, the vertical velocity of a large slope angle has a positive value in the end (t > 7 s). It is because the model box is finite in width (5H0), and soil particles falling to the wall begin to go upward about eight seconds later. Meanwhile, Figure 10 shows the cloud chart of vertical velocity from the SPH simulations of culvert-soil collapse with different slopes at t = 8.125 s, indicating that the positive velocity around the right boundary in Figure 8(c) is far larger than velocity in Figures 10(a) and 10(b).

Figure 11 compares the culvert-soil interactions with a far-away right boundary and those with bonded boundaries. It could be seen that the closed boundary restricts the movement of culvert, lateral/vertical velocity of flowing sand, and lateral interaction force greatly. It means that if the soil at both sides is nonliquefied or reinforced, the culverts suffer quite small lateral displacement and continue to function, which correspond with the common sense.

5. Conclusions

This paper proposes a solid-fluid coupled SPH model to study the culvert-soil system under large deformation. In this model the basic concepts of the SPH method, the contact algorithm and the constitutive equation for large deformation of liquefied sand are presented. Besides, the momentum equation is added with the fluid pressure of liquefied soil on culverts.

Then, the feasibility and accuracy of the model have been validated by dam-break tests and flowing sand tests, from which the flowing distance and impact force express a good agreement with experimental data.

After that, the presented SPH model is employed to analyze the lateral force and movement of culverts buried in the flowing sand, in which the flowing behavior of liquefied sand is also involved.

Following conclusions could be drawn from the calculated results:

(1) The lateral force of flowing sand on culvert attenuates rapidly at the beginning stage, and the time consumption from the maximum lateral force to minimum is not affected by the change of culvert size and buried depth;

(2) For a more inclined surface, the time of the appearance of maximum lateral force on the culvert becomes earlier, the run-out distances of the culvert, and the sand flowing velocity increase greatly;

(3) The closed right boundary has a rigid restriction to the lateral spread of sand profile and, in the flowing process, the relatively large and deep buried culverts suffer the larger maximum lateral force.

Data Availability

Previously reported chart data were used to support this study and are available at DOI: 10.1007/s10346-011-0285-5 and 10.1098/rsta.1952.0006. These prior studies (and datasets) are cited at relevant places within the text as references [26, 35].

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This work has been supported by the National Key Research and Development Program of China (Grant no. 2016YFC0800205), the National Key Basic Research Program of China (973 Program) (Grant no. 2015CB057901), the National Natural Science Foundation of China (NSFC Grant nos. 41630638 and 51808192), the Natural Science Foundation of Jiangsu Province (Grant no. BK20170887), and the Postdoctoral Science Foundation of China.