Study on Two-Phase Fluid-Solid Coupling Characteristics in Saturated Zone of Subgrade Considering the Effects of Fine Particles Migration

Abstract

The fluid seepage in saturated zone of subgrade promotes the migration of fine particles in the filler, resulting in the change of pore structure and morphology of the filler and the deformation of solid skeleton, which affects the fluid seepage characteristics. Repeatedly, the muddy interlayer, mud pumping, and other diseases are finally formed. Based on the theory of two-phase seepage, the theory of porous media seepage, and the principle of effective stress in porous media, a two-phase fluid-solid coupling mathematical model in saturated zone of subgrade considering the effects of fine particles migration is established. The mathematical model is numerically calculated with the software COMSOL Multiphysics^®. The two-phase seepage characteristics and the deformation characteristics of the solid skeleton in saturated zone of the subgrade are studied. The research results show that the volume fraction of fine particles first increases then decreases and finally becomes stable with the increase of time, due to the continuous erosion and migration of fine particles in saturated zone of the subgrade. The volume fraction of fine particles for the upper part of the subgrade is larger than that for the lower part of the subgrade. The porosity, the velocity of fluid, the velocity of fine particles, and the permeability show a trend of increasing first and then stabilizing with time; the pore water pressure has no significant changes with time. The vertical displacement increases first and then decreases slightly with the increase of time, and finally tends to be stable. For the filler with a larger initial volume fraction of fine particles, the maximum value of the volume fraction of fine particles caused by fluid seepage is larger, and the time required to reach the maximum value is shorter. It can be concluded that the volume fraction of fine particles in the subgrade filler should be minimized on the premise that the filler gradation meets the requirements of the specification in actual engineering.

1. Introduction

Filler for railway subgrade is mostly coarse-grained soil materials, mixed with a certain proportion of fine-grained soil, and has complex pore structure characteristics. When atmospheric precipitation, surface water, or snow melt infiltrate into subgrade and fail to be discharged in time, the local area of the subgrade will be saturated. Under the action of external loads, excess pore water pressure will be generated in saturated zone of the subgrade, and the dissipation of excess pore water pressure will form a high-impact transient fluid, which will drive fine particles to seepage in the pores of the filler [1,2,3]. In the process of fluid seepage, there are complex effects of erosion and migration of fine particles, which will lead to the change of pore structure and morphology of the filler. At the same time, the action of external loads and the change of the pore structure of the filler will cause deformation of the solid skeleton, which will further change the fluid seepage channel and affect the fluid seepage characteristics, and ultimately induce a series of diseases such as mud interlayer, mud pumping, and so on [4,5]. It is of great significance to deeply explore the fluid seepage characteristics and solid skeleton deformation characteristics considering the effects of fine particles migration and the change of hydraulic characteristics of coarse-grained soil to analyze the disaster-causing mechanism of the above diseases.

At present, there are many studies on the model of fine particles migration and the mechanism of fine particles migration both theoretically and experimentally, and have achieved fruitful results. J. F. Schijven [6] and R. Kretzschmar [7] have established a classical model for fine particles migration including the conservation equation of fine particles quantity and the kinetic equation of deposition. S. A. Bradford et al. [8] proposed a fine particles migration model considering the effects of sedimentation and sieving. On the basis of the traditional deep filtration model, J. Zhao et al. [9] applied the law of mass conservation, and established a particles migration model that considers the influences factors of seepage reduction coefficient, the velocity reduction coefficient, the change of porosity, and the suspended particles capture probability. The research results of I. Alobaidi [4,5] showed that the saturated interlayer water inside subgrade produces excess pore water pressure under the action of cyclic loading, which causes the vertical seepage of the water and the liquefied fine particles. B. Kermani [10] believed that the excess pore water pressure generated inside subgrade under cyclic loading drives the migration of fine particles, and the migration velocity of fine particles was analyzed. D. Liu, J. P. Herzig, and other scholars [11,12,13,14,15] studied the influences of fine particles concentration and particles size on the characteristics of fine particles migration in porous media. Q. S. Liu [16] conducted an experimental study on the migration-deposition characteristics of suspended particles with different sizes in porous media, and divided the migration-deposition types of suspended particles into three types including filter cake filter type, migration-deposition type, and free movement type.

On the other hand, the current researches on fluid-solid coupling characteristics have achieved rich research results. The coupling theory of water and soil (rock) was first used to solve the consolidation problem of soil. Terzaghi [17] studied the coupling phenomenon of fluid-solid deformation, and first proposed the concept of effective stress and established a one-dimensional consolidation model of soil. Subsequently, Biot [18,19] established a relatively complete three-dimensional consolidation theory. Since the establishment of Biot’s three-dimensional consolidation theory, researchers have introduced different material constitutive relations and applied them to a wider range of fields. Zienkiewicz [20] considered nonlinear materials and proposed the generalized Biot formula on the basis of Biot’s three-dimensional consolidation theory. In 1982, Savage and Bradok [21] extended the Biot’s three-dimensional consolidation theory to anisotropic elastic porous media. Zimmer-manl [22] introduced the dimensionless parameters to express the coupling effect of seepage and stress, and further extended the theoretical model of the coupling of multiphase saturated seepage in porous media. Tahmasebi [23,24,25,26,27,28] studied the effects of relative permeability, initial water saturation, residual oil saturation, and other factors on the interaction of multiphase fluid and solids through different new methods such as the volume of fluid (VOF) method, the coupled method of discrete element method (DEM) and volume of Fluid (VOF), and the coupled method of discrete element method (DEM) and computational fluid dynamics (CFD). These research results have important guiding significance for modeling the interactions between the particles and fluid in a multiphase fluid system. In 1977, Z. J. Shen [29] first combined Biot’s consolidation theory with finite element method and applied it to slope stability analysis. Y. Liang [30] obtained a fluid-solid coupling model of piping failure based on the mass conservation equation and the continuity equation of seepage. Z. Zhou [31] conducted artificial rainfall simulation tests on soil-rock mixture slopes, and studied the permeability of soil-rock mixture. On the basis of fluid mechanics and soil mechanics, W. Chen [32] derived the equilibrium equation of porous media and the continuity equation of pore fluid considering the coupling of stress field and seepage field, and established a fluid-solid coupling model.

However, the above researches on the characteristics of fine particles migration and fluid-solid coupling are mostly separate, and seldom involve the influence of fine particles migration on the hydraulic characteristics of coarse-grained soils. There is a lack of systematic research on the coupling effect between fine particles migration, fluid seepage, and deformation characteristics of solid skeleton.

Therefore, based on the theory of two-phase seepage, the theory of porous media seepage, and the principle of effective stress in porous media, this paper establishes a two-phase fluid-solid coupling mathematical model considering the effects of fine particles migration. The mathematical model is numeralized by the software COMSOL Multiphysics^®. The characteristics of fine particles migration, the two-phase seepage and the deformation of solid skeleton are studied.

2. Two-Phase Fluid-Solid Coupling Mathematical Model Considering the Effects of Fine Particles Migration

2.1. Basic Assumption

Filler for railway subgrade is mostly coarse-grained soil materials, mixed with a certain proportion of fine-grained soil, and has complex pore structure characteristics. Therefore, the saturated coarse-grained filler can be regarded as a saturated porous medium composed of solid medium (coarse particles and fine particles), fluid medium in pores, and fine particles that erode into the fluid medium that can migrate with the fluid (hereinafter referred to as fine particles). The structure of saturated porous media is shown in Figure 1. In order to establish a two-phase fluid-solid coupling mathematical model in saturated zone of subgrade considering the effects of fine particles migration, the following assumptions are made:

(1)

The coarse-grained filler is uniform and isotropic linear elastomer, and satisfies the assumption of small strain;

(2)

The seepage of water in the pores is caused by the pore water pressure gradient and follows Darcy’s law;

(3)

The influence of temperature and chemical action on the properties of subgrade filler is not considered;

(4)

Solid particles and water are incompressible;

(5)

The movable fine particles are spheres with the same radius (d_p), and the influence of fine particles on the fluid seepage characteristics and the force between fine particles are ignored.

2.2. Establishment of Two-Phase Fluid-Solid Coupling Mathematical Model Considering the Effects of Fine Particles Migration

2.2.1. Equations of Seepage Field

(1) Basic definitions

As shown in Figure 1, the micro-model for coarse-grained soil is composed of three parts, which includes solid medium ( $s$, $d{m}_s$ is mass of the solid and $d{V}_s$ is volume of the solid), fluid medium ($l$, $d{m}_l$ is mass of the solid and $d{V}_l$ is volume of the solid), and fine particles ($p$, $d{m}_p$ is mass of the solid and $d{V}_p$ is volume of the solid).

Porosity of porous media is defined as:

$$\varphi =\frac{d{V}_V}{dV}=\frac{d{V}_l+d{V}_p}{dV}$$

(1)

where: $\varphi$ is the porosity of the micro-model; $d{V}_V$ is the pore volume of the micro-model (${\mathrm{m}}^3$); $dV$ is the total volume of the micro-model (${\mathrm{m}}^3$).

The volume fraction of fine particles is:

$$C=\frac{d{V}_p}{d{V}_V}$$

(2)

where: $C$ is the volume fraction of fine particles.

(2) Equations of mass conservation

According to the law of mass conservation of fluid, when considering the two-phase (fluid and fine particles) seepage in porous media, the mass conservation equations of fluid phase and fine particles phase can be expressed as:

$$\frac{\partial }{\partial t}[\varphi (1-C){\rho }_l]+\nabla [(1-C){\rho }_l{v}_{ls}]=q_l$$

(3)

$$\frac{\partial }{\partial t}(\varphi C{\rho }_p)+\nabla (C{\rho }_p{v}_{ps})=\dot{m}$$

(4)

where: ${\rho }_l$ is the inherent density of the fluid ($\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$); ${\rho }_p$ is the inherent density of fine particles ($\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$); $v_{ls}$ is the velocity vector of the fluid relative to the solid ($\mathrm{m}/\mathrm{s}$); $q_l$ is the source (sink) item of the fluid ($\mathrm{k}\mathrm{g}/{\mathrm{m}}^3\mathrm{s}$); $v_{ps}$ is the velocity vector of fine particles relative to solid ($\mathrm{m}/\mathrm{s}$); $\dot{m}$ is the mass rate of fine particles migration ($\mathrm{k}\mathrm{g}/{\mathrm{m}}^3\mathrm{s}$); $\nabla =\frac{\partial }{\partial x}i+\frac{\partial }{\partial y}j+\frac{\partial }{\partial z}k$ is Hamiltonian operator.

(3) Equations of motion

The seepage of fluid conforms to Darcy’s law, as follows:

$$v_{ls}=v_l-v_s=-\frac{k}{{\mu }_l}(\nabla p_l-{\rho }_{l}g)$$

(5)

where: $v_l$ is the velocity vector of the fluid ($\mathrm{m}/\mathrm{s}$); $v_s$ is the velocity vector of solid ($\mathrm{m}/\mathrm{s}$); $k$ is the permeability of porous media (${\mathrm{m}}^2$); $p_l$ is the pore pressure of the fluid ($\mathrm{P}\mathrm{a}$); ${\mu }_l$ is the dynamic viscosity coefficient of the fluid ($\mathrm{P}\mathrm{a}\cdot \mathrm{s}$); $g$ is the gravitational acceleration vector ($\mathrm{N}/\mathrm{k}\mathrm{g}$).

The force of fine particles in liquid is complex. The fine particles is mainly affected by effective gravity, drag force, additional mass force, pressure gradient force, Basset force, Magnus force, Saffman force, and other force, among which additional mass force and pressure gradient force, Basset force, Magnus force, Saffman force have relatively small influence on the movement of fine particles and can be ignored. Therefore, considering the effect of effective gravity and drag force on fine particles, the schematic diagram of the force on a single fine particle is shown in Figure 2.

We use $v_r$ to represent the difference between the velocity vector of the fluid and the velocity vector of the fine particles ($v_p$). Assuming that the velocity of the fine particles differs little from that of fluid in the horizontal direction, and the velocity difference in the vertical direction can be expressed as:

$$v_{rz}=v_{lz}-v_{pz}$$

(6)

When a single particle is in force equilibrium in the fluid, there are:

$$({\rho }_p-{\rho }_l)\frac{\pi d_p^3}{6}g=C_D\frac{\pi d_p^2}{4}\frac{{\rho }_l{v_{rz}}^2}{2}$$

(7)

then:

$$v_{rz}=\sqrt{\frac{4({\rho }_p-{\rho }_l)d_{p}g}{3C_D{\rho }_l}}$$

(8)

where, $d_p$ is the diameter of fine particles ($\mathrm{m}$);$C_D$ is the coefficient of drag force, which is related to the Reynolds number of fine particles.

$$\mathrm{Re}=\frac{d_p{\rho }_l{u}_l}{{\mu }_l}$$

(9)

$$\{\begin{array}{l}{C}_D=\frac{24}{\mathrm{Re}},\mathrm{Re}<2, \mathrm{Laminar} \mathrm{flow} \mathrm{zone};\\ C_D=\frac{18.5}{{\mathrm{Re}}^{0.6}},2<\mathrm{Re}<500, \mathrm{Transition} \mathrm{zone};\\ C_D\approx 0.44,500<\mathrm{Re}<2\times {10}^5, \mathrm{Turbulence} \mathrm{zone};\\ C_D\approx 0.11,\mathrm{Re}>2\times {10}^5, \mathrm{Turbulent} \mathrm{boundary} \mathrm{zone};\end{array}$$

(10)

According to Equations (6)~(10), the motion equations of fine particles can be obtained:

$$\{\begin{array}{l}{v}_{psx}=-\frac{k}{{\mu }_l}\frac{\partial p_l}{\partial x}\\ v_{psy}=-\frac{k}{{\mu }_l}\frac{\partial p_l}{\partial y}\\ v_{psz}=-\frac{k}{{\mu }_l}(\frac{\partial p_l}{\partial z}-{\rho }_{l}g)-\sqrt{\frac{4({\rho }_p-{\rho }_l)d_{p}g}{3C_D{\rho }_l}}\end{array}$$

(11)

2.2.2. Equations of Stress Field

According to the principle of effective stress in saturated porous medium, there are:

$${\sigma }_{ij}={\sigma }_{ij}^{\prime }+\varphi p_l{\delta }_{ij}$$

(12)

where, ${\sigma }_{ij}$ is the total stress tensor; ${{\sigma }^{\prime }}_{ij}$ is the effective stress tensor; ${\delta }_{ij}$ is the Kronecker symbol,

$${\delta }_{ij}=\{\begin{array}{l}1,i=j\\ 0,i\ne j\end{array}.$$

According to the assumption that the solid skeleton is uniform and isotropic linear elastomer, and satisfies the assumption of small strain, its constitutive equation is:

$${{\sigma }^{\prime }}_{ij}=\lambda {\epsilon }_V{\delta }_{ij}+2\mu {\epsilon }_{ij}$$

(13)

where, ${\epsilon }_V$ is the volumetric strain of the solid skeleton, ${\epsilon }_V={\epsilon }_x+{\epsilon }_y+{\epsilon }_z$, and ${\epsilon }_x,{\epsilon }_y,{\epsilon }_z$ is the axial strain in the three directions of $x$, $y$, and $z$, respectively; ${\epsilon }_{ij}$ is the strain of the solid skeleton; $\lambda , \mu$ is the constant of Lamé,

$$\lambda =\frac{Ev}{(1+v)(1-2v)}$$

(14)

$$\mu =\frac{E}{2(1+v)}$$

(15)

where, $E$ is Young’s elastic modulus; $v$ is Poisson’s ratio.

The displacement of the solid skeleton is expressed in terms of $u_s$, and the relation between strain and displacement is:

$${\epsilon }_{ij}=\frac{1}{2}(u_{sj,i}+u_{si,j})$$

(16)

Then, the stress balance equation of porous media can be expressed as:

$${\sigma }_{ij,j}+R_i=0$$

(17)

where, $R_i$ is the volume force of the porous medium ($\mathrm{N}$);

$$R_i=\rho g=\left[\right(1-\varphi ){\rho }_{\mathrm{s}}+\varphi (1-C){\rho }_l+\varphi C{\rho }_p]g$$

(18)

where, $\rho$ is the average density of the porous medium ($\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$); ${\rho }_s$ is the inherent density of the solid ($\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$).

Then, the governing equation of the stress field of porous media can be obtained from Equations (12), (13), and (16)~(18).

$$(\lambda +\mu )\nabla {\epsilon }_V+\mu {\nabla }^2{u}_{si}+\nabla (\varphi p_l{\delta }_{ij})+R_i=0$$

(19)

2.2.3. Migration Equation of Fine Particles

According to the references [33,34,35,36], considering the limit value of the migration of fine particles [7,37,38,39]. The migration equation of fine particles is as follows:

$$\dot{m}={\rho }_p\zeta ({\varphi }_{\max }-\varphi )(C-\frac{C^2}{C_{cr}})\left|v_l\right|$$

(20)

where: ${\varphi }_{\max }$ is the maximum porosity, which refers to the limit value of porosity change in porous media caused by the migration of fine particles; $C_{cr}$ is the critical volume fraction of fine particles when the migration reaches equilibrium; $\left|v_l\right|$ is the modulus of the velocity vector of the fluid ($\mathrm{m}/\mathrm{s}$); $\zeta$ is the erosion coefficient of fine particles (${\mathrm{m}}^{-1}$).

2.2.4. Auxiliary Equation

(1) Porosity evolution equation

The change in pore volume is mainly caused by the erosion and migration of fine particles and the deformation of solid skeleton. For the erosion and migration of fine particles, the rate of change in pore volume is equal to the volume of eroded fine particles per unit time, then:

$$\frac{\partial (d{V}_V)}{\partial t}=\frac{\dot{m}dV}{{\rho }_p}$$

(21)

Which can be simplified as:

$$\frac{\partial \varphi }{\partial t}=\frac{\dot{m}}{{\rho }_p}$$

(22)

For pore volume change caused by deformation of solid skeleton, the porosity equation proposed in reference [40] is as follows:

$$\varphi =\frac{{\varphi }_0+{\epsilon }_V}{1+{\epsilon }_V}$$

(23)

where, ${\varphi }_0$ is the initial porosity.

Equation (23) can be written as:

$$\frac{\partial \varphi }{\partial t}=\frac{1-{\varphi }_0}{{(1+{\epsilon }_V)}^2}\frac{\partial {\epsilon }_V}{\partial t}$$

(24)

Then, based on the erosion and migration of fine particles and the deformation of solid skeleton, the porosity evolution equation can be expressed as:

$$\frac{\partial \varphi }{\partial t}=\frac{\dot{m}}{{\rho }_p}+\frac{1-{\varphi }_0}{{(1+{\epsilon }_V)}^2}\frac{\partial {\epsilon }_V}{\partial t}$$

(25)

(2) Relationship between permeability and porosity

The Kozeny-Carman equation was used to describe the relationship between permeability and porosity [41,42,43]:

$$k=k_0{(\frac{\varphi }{{\varphi }_0})}^3{(\frac{1-{\varphi }_0}{1-\varphi })}^2$$

(26)

where: $k_0$ is the initial permeability (${\mathrm{m}}^2$).

Equations (3)~(5), (11), (19), (20), (25) and (26), constitute the two-phase fluid-solid coupling mathematical model considering the effects of fine particles migration.

2.2.5. Boundary Conditions and Initial Conditions

The coupling mathematical model mentioned above is highly nonlinear, boundary conditions and initial conditions need to be given in the process of solution. The pore water pressure, the volume fraction of fine particles, and the displacement of solid skeleton are taken as basic variables.

Boundary conditions: the boundary conditions usually have the first type, also known as Dirichlet boundary; and the second type, also known as Neumann boundary.

The first boundary conditions can be expressed as:

$p_l={\overline{p}}_l$ on ${\Gamma }_{p_l}$;

$C=\overline{C}$ on ${\Gamma }_C$;

$u_s={\overline{u}}_s$ on ${\Gamma }_{u_s}$.

The second boundary conditions can be expressed as:

$q_l=-v_{ls}\cdot n$ on ${\Gamma }_{q_l}$;

$q_C=-C{v}_{ps}\cdot n$ on ${\Gamma }_C$;

${\sigma }_{ij}\cdot n=\overline{t}$ on ${\Gamma }_{{\sigma }_{ij}}$.

Initial conditions:

$p_l=p_l^0$, $C=C^0$, $u_s=u_s^0$ at t = 0.

2.2.6. The Solution of the Mathematical Model

The above mathematical model was solved by using the coupling solid mechanics module and porous medium multiphase flow module in the software COMSOL Multiphysics^®. Taking the volume fraction of fine particles ($C$), the pore pressure of the fluid ($p_l$), and the displacement of the solid skeleton ($u_s$) as the basic unknowns, the other unknowns such as the velocity vector of the fluid $v_l$, the velocity vector of the fine particles ($v_p$), the mass rate of fine particles migration ($\dot{m}$), the porosity ($\varphi$), and the permeability ($k_0$)can be calculated by the basic unknowns.

3. Validation of the Coupling Model

In this section, the software COMSOL Multiphysics^® is used to solve the mathematical model. Moreover, the solution results were compared with those in reference [33] to verify the accuracy of the coupling model in this paper.

3.1. The Establishment of Numerical Model

Reference [33] studied the problem of sand production in the wellbore during oil extraction. A two-dimensional calculation model is adopted. The radius of the wellbore (inner boundary) is 0.1 m, and the outer boundary is 5 m. The calculation model is shown in Figure 3.

The standard atmospheric pressure value is taken as the calculation reference value. The initial porosity of the model was 0.25, and the initial volume fraction of fine particles was 0.001. The boundary conditions adopted by the model are as follows: the fluid pressure at inner boundary is 5 MPa, and the fluid pressure at outer boundary is 8 MPa. The inner boundary is the outflow boundary for fine particles, and the external stress acting on the outer boundary is 20 MPa. The physical parameters are shown in

Table 1. Physical parameters for numerical calculation.

Parameter	Value
Wellbore(inner boundary) radius $r_a$/(m)	0.1
Outer boundary radius $r_o$/(m)	5.0
Initial porosity ${\varphi }_0$	0.25
Maximum porosity ${\varphi }_{\max }$	1
Initial volume fraction of fine particles $C_0$	1.0 × 10−3
Critical volume fraction of fine particles $C_{cr}$	0.30
Initial permeability $k_0$/(m2)	3.73 × 10−13
The dynamic viscosity coefficient of the fluid ${\mu }_l$/(Pa·s)	4.2 × 10−3
The density of fluid ${\rho }_l$/(kg/m3)	1000
The density of fine particles ${\rho }_p$/(kg/m3)	2650
The erosion coefficient of fine particles $\zeta$/(m−1)	5.0
Fluid pressure at inner boundary $p_{l-a}$/(MPa)	5.0
Fluid pressure at outer boundary $p_{l-o}$/(MPa)	8.0
Elastic model $E$/(MPa)	2.0 × 103
Poisson’s ratio $\nu$	0.3
External stress ${\sigma }_0$/(MPa)	20

.

3.2. Verification and Analysis of Calculation Results

The evolution curves of volume fraction of fine particles, porosity, velocity of fluid, and radial displacement at the inner boundary calculated according to the mathematical model in this paper are shown in Figure 4, Figure 5, Figure 6 and Figure 7 respectively, and the calculation results in this paper are compared with those in reference [33]. A comprehensive analysis of Figure 4, Figure 5, Figure 6 and Figure 7 shows that the evolution curves of volume fraction of fine particles, porosity, velocity of fluid, and radial displacement calculated by the mathematical model established in this paper are consistent to the results in reference [33]. Because of the difference between the porosity evolution equation, the permeability equation in this paper, and those in reference [33], there is a certain error between the calculation results and those in reference [33]. However, the error is within the acceptable range. This also verifies the correctness of the model in this paper.

4. Analysis of Fluid-Solid Coupling Characteristics of Saturated Coarse-Grained Filler Considering the Effects of Fine Particles Migration

4.1. The Establishment of Numerical Model

According to “Code for Design of Heavy Haul Railway (TB10625-2017)” [44] in China and referring to a typical section structure of the Shuo-huang Railway subgrade, a two-dimensional calculation model is established as shown in Figure 8. The structure of the subgrade from top to bottom is as follows: surface layer of the subgrade (the thickness is 0.6 m), bottom layer of the subgrade (the thickness is 1.9 m). The slope gradient of the subgrade is 1:1.5. The model sets four observation points, as shown in Figure 8. According to references [45,46,47,48], we use the free mesh triangular elements and an ultra-fine mesh (the number of mesh is 922, the minimum mass for the mesh is 0.55 and the average mass for the mesh is 0.89) is set in this research. The mesh of the calculation model is shown in Figure 9a, and the influence of mesh number on the velocity of fluid and the influence of mesh number on the vertical displacement is shown in Figure 9b. The results indicate that, the calculation results of models with the mesh number of 138, 303, and 922 are different. When the mesh number is increased to 922, the calculation results are little different from those of model with the mesh number 3566. Therefore, the ultra-fine mesh in this paper is accepted.

4.2. Mechanical Parameters of Subgrade Filler

The filler of the subgrade is made of graded crushed stone. In order to analyze the influence of volume fraction of fine particles on the fluid-solid two-phase coupling characteristics of the filler, three different volume fractions of fine particles are selected, and the gradation of the graded crushed stone filler meet the requirements of the specification [49]. The physical parameters are shown in

Table 2. Physical parameters for numerical calculation of subgrade model.

Parameter	Value
Initial porosity ${\varphi }_0$	0.33
Maximum porosity ${\varphi }_{\max }$	0.50
Initial volume fraction of fine particles $C_0$	0.01/0.03/0.07
Critical volume fraction of fine particles $C_{cr}$	0.30
Initial permeability $k_0$/(m2)	7.80 × 10−12
The dynamic viscosity coefficient of the fluid ${\mu }_l$/(Pa·s)	1.0 × 10−3
The density of fluid ${\rho }_l$/(kg/m3)	1000
The density of fine particles ${\rho }_p$/(kg/m3)	2380
The erosion coefficient of fine particles $\zeta$/(m−1)	1.0
Fluid pressure $p_{l-in}$/(kPa)	40.0
Elastic model $E$/(MPa)	120
Poisson’s ratio $\nu$	0.3

.

4.3. Boundary Conditions and Initial Conditions

The initial conditions for the calculation model conclude that the pore water pressure is the hydrostatic pressure; the initial volume fractions of fine particles are 0.01, 0.03, and 0.07, respectively; the initial porosity of the subgrade is 0.3; and the initial displacement is 0.

It is assumed that the subgrade in the calculated area is saturated, and the dissipation path of pore water pressure is from bottom to top. The seepage boundary conditions for the calculation model are as follows: the lower boundary AD is the boundary of fluid pressure, and the fluid pressure at the boundary AD is set as 40 kPa according to the reference [1,2,3]. The fluid pressure at upper boundary BC is 0 kPa. The slopes are usually protected, so AB and CD are assumed to be impermeable boundaries. The boundary conditions corresponding to the migration of fine particles are as follows: the boundary AD is the inflow boundary for fine particles carried by the fluid, and the volume fractions of fine particles at the boundary AD are 0.01, 0.03, and 0.07, respectively. The upper boundary BC is the outflow boundary of fine particles. AB and CD are no flux boundary. The displacement boundary conditions are as follows: the boundary AD is set as a fixed boundary, the remaining boundaries are free boundaries.

4.4. Calculation Results and Analysis

4.4.1. The Variation of the Volume Fractions of Fine Particles

Figure 10 shows the evolution curves of the volume fractions of fine particles in saturated zone of railway subgrade. We can see from Figure 10 that:

(1)

With the increase of time, the volume fractions of fine particles at the observation points 1 to 3 first increased and then decreased: at initial moment, the volume fraction of fine particles was 0.01; when time is between 1100 min and 1400 min, the volume fraction of fine particles at different observation points in the subgrade reaches the maximum. Subsequently, the volume fractions of fine particles decreases gradually to about 0.01. This indicates that the fine particles are eroded by the fluid and migrated with the fluid, which leads to the increase of the volume fractions of fine particles in the subgrade. As time goes by, fine particles are continuously eroded and flow out, leaving skeleton particles and larger particles in the subgrade which are difficult to migrate. Therefore, the volume fractions of fine particles decrease gradually.

(2)

For observation points 1 to 3, the volume fraction of fine particles at observation point 3 first reaches the maximum value, and then the volume fraction of fine particles at observation point 2 and 1 successively reaches the maximum value. The maximum volume fraction of the fine particles observed at point 1 is greater than that observed at point 2, the maximum volume fraction of the fine particles observed at point 2 is greater than that observed at point 3. The reasons are as follows: the seepage path of the fluid in the subgrade is from bottom to top, and the fluid erodes fine particles, then carries fine particles from the lower part of the subgrade to the upper part of the subgrade. Therefore, the volume fraction of fine particles at the observation point 3 first reaches its maximum value. However, the fluid continuously erodes fine particles and the eroded fine particles migrate to the upper part of the subgrade and flow out of the subgrade from the upper boundary. Therefore, the maximum volume fraction of fine particles at observation point 1 is larger than that at other observation points.

(3)

As the observation point 4 is located at the lower boundary, it is the inflow boundary of fine particles, so the volume fraction of fine particles at the observation point 4 does not change.

4.4.2. The Variation of Porosity

Figure 11 shows the evolution curves of porosity in saturated zone of railway subgrade, in which subfigure (b) is the detailed curves for the curves in the box of subfigure (a). It can be seen from Figure 11 that:

(1)

The variation trend of porosity at the four observation points is same. First, it increases slowly, then it increases rapidly to a certain value (about 0.5, i.e., the maximum porosity), and then it is basically stable.

(2)

For observation points 1 and 2, the time for the porosity to reach the maximum porosity is same; for observation point 3, the time for the porosity to reach the maximum porosity is slightly later than that of observation points 1 and 2; while for observation point 4, it has the latest time to reach the maximum porosity.

4.4.3. The Variation of Pore Water Pressure

Figure 12 shows the distribution cloud map of pore water pressure in saturated zone of railway subgrade, and Figure 13 shows the evolution curve of pore water pressure at the observation points 3. From these two pictures we can see that:

(1)

The pore water pressure at the lower boundary of the subgrade is 65 kPa, and the pore water pressure at the upper boundary of the subgrade is 0 kPa. The pore water pressure is evenly distributed from the bottom of the subgrade to the top of the subgrade.

(2)

With the increase of time, the pore water pressure at point 3 shows a trend of first decreasing, then increasing, and then decreasing, and finally remains stable, but the range of its increase or decrease is small. This reason is that with the erosion and migration of fine particles in the subgrade, the porosity and permeability of each point inside the subgrade are different (see Figure 10 and Figure 14 for details), thus affecting the distribution of pore water pressure. Therefore, the pore water pressure at the same point is not constant.

4.4.4. The Variation of the Velocity of Fluid and the Velocity of Fine Particles

Figure 14 shows the evolution curves of the velocity of fluid and the velocity of fine particles in saturated zone of railway subgrade. From Figure 14 we can see:

(1)

The variation trend of fluid velocity is basically the same as that of fine particles migration velocity, showing a trend of basically unchanged at first, then rapidly increasing, and finally tending to be stable.

(2)

The velocity of fluid is always greater than the velocity of fine particles, and the velocity difference between these two is becoming more and more obvious with the increase of fluid velocity. The reason is that with the increase of the fluid velocity, the drag coefficient decreases, and the drag force of the fine particles is reduced, which hinders the movement of the fine particles to some extent.

Figure 15 shows the relation curves between fluid velocity and porosity at different observation points. It can be seen from Figure 15 that the difference between the fluid velocity and porosity at different observation points is small. The fluid velocity increases with the increase of porosity, and shows an exponential change.

4.4.5. The Variation of Permeability

Figure 16 shows the evolution curves of permeability in saturated zone of railway subgrade. It can be seen from Figure 16 that the permeability increases rapidly to the maximum after a relatively slow growth, and finally remains stable. For observation points 1~3, the maximum permeability is about 4.9 × 10⁻¹¹ m², which is an increase of about 6.3 times higher than the initial permeability.

4.4.6. Effects of Initial Volume Fraction of Fine Particles on the Evolution Curves of the Volume Fraction of Fine Particles

Figure 17 shows the evolution curves of the volume fractions of fine particles with different initial fine particles volume fraction at point 2. From Figure 17, we know that the larger the initial fine particles volume fraction is, the larger the maximum value of fine particles volume fraction is, and the shorter the time needed to reach the maximum value. Combined with Equation (20), the larger the initial fine particles volume fraction is, the larger the mass rate of fine particles migration will be; therefore, the shorter the time needed to reach the maximum fine particles volume fractions will be. It can be concluded that the fine particles content in the filler of subgrade should be strictly controlled, and the fine particles content should be reduced as far as possible on the premise that the filler grade meets the requirements of the specification in practical engineering.

4.4.7. The Variation of Effective Stress

Figure 18 shows the contours of vertical effective stress at the steady state of fluid seepage. Figure 19 shows the evolution curves of vertical effective stress at different observation points in saturated zone of railway subgrade. It can be seen from Figure 18 that the vertical effective stress is distributed symmetrically on the central axis of the subgrade. With the increase of depth, the vertical effective stress decreases. It can be seen from Figure 19 that vertical effective stress at points 2 and points 3 increases first and then decreases slightly with the increase of time, and finally tends to be stable. In combination with Figure 13, the pore water pressure increases at a time of about 1000–2000 min, so the effective stress decreases.

4.4.8. The Variation of the Vertical Displacement

Figure 20 shows the contours of vertical displacement at the steady state of fluid seepage. We can see that the vertical displacement is distributed symmetrically on the central axis of the subgrade. At the same depth, the vertical displacement at the central axis of the subgrade is slightly larger than the vertical displacement at the edge of the subgrade.

Figure 21 shows the evolution curves of vertical displacement at different observation points in saturated zone of railway subgrade. It can be seen from Figure 21 that the vertical displacement at different observation points increases first and then decreases slightly with the increase of time, and finally tends to be stable. The vertical displacement of the upper part of the subgrade is greater than that of the lower part of the subgrade. The final vertical displacement at observation point 1 is about 1.4 mm, and the final displacement at observation point 3 is about 0.5 mm.

4.4.9. The Relationship between Porosity Change and Accumulative Volume Fraction of Fine Particles/Displacement

Figure 22a shows the relation curves of porosity change and accumulative volume fraction of fine particles and (b) shows the relation curves of porosity change and vertical displacement. The positive value of porosity change represents the increase of porosity, and the negative value of porosity change represents the decrease of porosity. It can be seen from Figure 22 that the accumulative volume fraction of fine particles has a significant impact on the change of porosity: with the increase of the accumulation volume fraction of fine particles, porosity change increases gradually, when the accumulation volume fraction of fine particles reaches about 10, the porosity change reaches a plateau value. The vertical displacement has little influence on the change of porosity. With the change of vertical displacement, the porosity first decreases slightly and then increases slightly, and finally remains stable.

5. Conclusions

Based on the theory of two-phase seepage, the theory of porous media seepage, and the principle of effective stress in porous media, a two-phase fluid-solid coupling mathematical model in saturated zone of subgrade considering the effects of fine particles migration is established. The mathematical model is numerically calculated with the software COMSOL Multiphysics^®. The two-phase seepage characteristics in saturated zone of the subgrade and the deformation characteristics of the solid skeleton are studied, and the following conclusions are obtained:

(1)

Under the action of fluid seepage, the erosion and migration of the fine particles in saturated zone of subgrade gradually occurred, and the volume fraction of fine particles increased with the increase of time, and as the time goes on, only skeleton particles and particles with large sizes are left in the subgrade, the volume fraction of fine particles gradually decreased to the initial value. The constant erosion of fine particles inside the subgrade and the constant migration of fine particles with fluid from the lower part of the subgrade to the upper part of the subgrade caused that the volume fraction of fine particles in the upper part of the subgrade is larger than that in the lower part of the subgrade at the same time.

(2)

With the increase of time, the porosity, the velocity of fluid, the velocity of fine particles, and the permeability of the subgrade show the trend of first increasing and then stabilizing; the pore water pressure in the subgrade has no significant changes with time.

(3)

The vertical displacement is distributed symmetrically on the central axis of the subgrade. With the increase of time, the vertical displacement increases first and then decreases slightly and finally tends to be stable.

(4)

For the filler with a larger initial volume fraction of fine particles, the maximum value of the volume fraction of fine particles in the subgrade caused by fluid seepage is larger, and the time required to reach the maximum value is shorter. It can be concluded that in actual engineering, the volume fraction of fine particles in the subgrade filler should be minimized on the premise so that the filler gradation meets the requirements of the specification.

The main purpose of this research is to explore the characteristics of the two-phase seepage and the deformation of solid skeleton in saturated zone of the subgrade considering the effects of fine particles migration, in order to explain the internal mechanisms of the mud pumping disease. Therefore, the dynamic characteristics of the subgrade filler and the fluid-solid two-phase seepage characteristics considering the effects of fine particles migration and sedimentation under train loads have more practical engineering significance, which will also be the focus of the next research.