Plane shock wave propagation in partially saturated soil

1 Introduction

In recent years, the scale of cities in various regions has continued to expand, and underground transportation networks such as subways and underground tunnels have grown. The distribution and mechanical characteristics of soils in different regions are usually not the same, so the geological environment that underground tunnels pass through is quite complicated, especially for underwater shield tunnels. As shown in Table 1, underwater shield tunnels often pass through a variety of soil layers (Hong, 2007;27:1–3+18.). For example, the Shanghai Yangtze River Tunnel has an average burial depth of 9 m, and the shallowest burial depth is only 6.8 m. It passes through grey silt clay, grey clay, and other soil layers (Yuan et al., 2008). The air quantity in silt clay is usually very small, and it varies with depth within a certain range (Jiang and Luo, 2007).

The peak stress generally decreases in the direction of wave propagation in soils, and the peak stress on the rigid wall is not larger than double the incident stress. These phenomena have been theoretically and experimentally validated on many occasions for unsaturated soils and saturated soils (Drake and Little, 1983; Department of the Army. Fundamentals of protective design for conventional weapons. Technical manual TM5-855-1. Washington, DC: Department of the Army;, 1986; Lyakhov, 1974; Seaman and Whitman, 1966). Nevertheless, for partially saturated soils in the natural state, unique phenomena can be observed under transient loading (Lyakhov, 1974; Lyakhov, 1964; Dontsov and Nakoryakov, 2001; Emery and Valent, 1995), and the content of air in soils is the most dominant factor affecting the behaviour; unfortunately, there is a lack of details related to such experiments.

Due to the increasing-hardening constitutive relation of partially saturated soil, the propagation law of shock waves in such soil is different from that in other geotechnical media. The wave velocity in partially saturated soil is usually higher than that in unsaturated soil under the same explosion load. The rise time of the shock wave in the near and middle areas of the explosion is very short due to the increasing-hardening constitutive relation, which causes the structure to easily experience strong vibration and even be destroyed (Mu et al., 2010). The interaction between the soil and structure may also lead to the accumulation of pore pressure near the structure, which would destroy the stability of the structure (Liang et al., 2018). Therefore, it is very meaningful to study the propagation law of shock waves in partially saturated soil.

Lyakhov (Lyakhov, 1964) experimentally obtained the relationship between the peak free-field stress and scaled range for explosions in sandy soils. The experiments were performed in fine-grained sands with different air contents, and the charges were buried at a sufficient depth of 2.5Q^1/3 to exclude the effect of the free surface. According to the viewpoint of G. M. Lyakhov, in water-saturated soils, stress waves always have the features of shock waves. The stress resulting from wave propagation does not depend on the orientation, and the soil behaves like a liquid. In unsaturated soil, the stress in the direction of propagation is higher than that in the normal direction. It must be emphasized that all maximum stresses in the experiment performed by Lyakhov were greater than 1.0 MPa for water-saturated sand. According to the experimental results (Lyakhov, 1964), in a water-saturated soil at normal incidence of the shock wave on a rigid obstacle, the stress increases 2 to 3.3 times depending on the kind of soil and the magnitude of the incident stress. For low incident stress, the increase is twofold, and as the incident stress increases, it may reach a value as high as the upper limit given above. Furthermore, the stress in the media is considerably influenced by the state of the air content in the soil.

Based on Lyakov’s work, Bloom (Bloom, 2006) studied the problem of shock propagation in the plastic gas zone and proposed that Lyakov’s results were very useful in research of shock waves in soils. Richart et al. (Richart et al., 1970) investigated the theory for the propagation of elastic compression waves in partially saturated soils and discussed the effects of small amounts of air on the elastic wave propagation velocity. Mesgouez et al. (Mesgouez et al., 2005) and Schanz and Cheng (Schanz and Cheng, 2000) studied the stress and particle displacement of elastic waves in porous media based on Biot (Biot, 1956) theory. Zhou et al. (Zhou et al., 2013) used numerical examples to investigate the propagation characteristics of elastic waves in nonuniform saturated soil, such as the wavenumber and dissipation. Li et al. (Li et al., 2018), Kim et al. (Kim et al., 2002; Kim et al., 2002) and Xiang-Wei and Xiang-Yue (Xiang-wei and Xiang-yue, 1989) used different theories to analyse wave propagation properties in saturated porous media, and the influence of parameters such as the saturation and permeability in saturated soil on the wave velocity was studied. Wang et al. (Wang et al., 2019) measured longitudinal wave velocity and a comprehensive in several sandstone samples with different degrees of saturation and different micro-pore structure characteristics by improved acoustic emission. The results show that the wave velociy has almost no relationship with the pore size in the medium, but with the increase of saturation, the wave speed increases obviously. The wave velocity is related to the wave impedance, so the increase of the saturation will directly lead to the increase of the wave impedance, thus affecting the characteristics of the shock wave in the medium. Wang (Wang, 2021) investigated the nonlinear dynamic response in saturated soil caused by shock wave propagation based on numerical and experimental results. Gang et al. (Gang et al., 2019) comparatively analysed the experimental data of four groups of partially saturated soils under plane impact loading. The results showed that the dynamic behaviour of the soil was influenced by the water content boundary conditions when the dry densities were the same. Yu et al. (Yu et al., 1982) and Zhao et al. (Zhao et al., 2004; Zhao et al., 2004) conducted experiments on the propagation of shock waves in silty saturated soil, saturated sand, and saturated clay in a certain area, and the inverse-attenuation phenomenon of the free-field stress in saturated soil was discovered. When the partially saturated soil was close to fully saturated, a slight change in the air content would have a substantial effect on the shock wave propagation characteristics in soil (Li et al., 2018; Zhao et al., 2007). Because the maximum deformation of the structure under incident loading is often determined by the peak stress, it is meaningful to calculate the peak stress around the tunnel by studying the relationship between peak stress and depth in partially saturated soil.

2 Constitutive model of partially saturated soil

2.1

2.1 P-α model of partially saturated soil

The law of shock wave propagation in a specific physical model can be intuitively and effectively obtained by a set of experiments. However, it is very difficult to precisely control the saturation of the sample, and the scale of the experiment is limited by the experimental scale, budget, etc. More importantly, the geological conditions vary substantially by region. These objective conditions make it difficult to propose the propagation law of shock waves in various geological conditions for experimental applications.

Bloom (Bloom, 2006) has proposed that Lyakov’s results were very useful in research of shock waves in soils. But the phenomenon that the gas content of partially saturated soil would decrease and the saturation will gradually increase under a shock load was not considered in the constitutive model. So Lyakhov's nonlinear elastic model was improved to a nonlinear model that considered volumetric plastic deformation in this paper.

The volume compression equations of air, water, and solid particles can be described by the Tait equation:

(1)

$$P_f=\frac{{\rho }_{i0}{C}_{i0}^2}{{\gamma }_i}\left[{\left(\frac{{\rho }_i}{{\rho }_{i0}}\right)}^{{\gamma }_1}-1\right]i=1,2,3$$

where P_f represents the hydrostatic pressure inside the medium, and the hydrostatic pressure acting on different components is constant. i = 1, 2, 3 represent air, water, and solid particles, respectively; C_i0 represents the wave velocity of the i-th component; and ρ_i0 and ρ_i represent the density before and after compression, respectively, of the i-th component. γ_i represents the isentropic index of the i-th component, where γ_1,2,3 = 1.4, 3.0, and 7.0.

The skeleton part is approximated by a bilinear stress–strain relationship, as shown in Fig. 1(a). The relationship of elasticity gives the following equation:

(2)

$$K=\frac{E\nu }{\left(1+\nu \right)\left(1-2\nu \right)}+\frac{2}{3}G=\frac{E}{3\left(1-2\nu \right)}$$

Close

Fig. 1

Sketch of the Compressive Stress-Strain Relationship.

Export to PPT

Where K represents the bulk modulus, E represents the elastic modulus, v represents the poisson's ratio, and G represents the shear modulus.

The volume compression relation of the skeleton in the elastic deformation stage can be easily obtained by Eq. (2).

After the explosion wave propagates, the volume residual deformation, which mainly comes from the dissolution of gas and the restacking of the skeleton structure, in partially saturated soil is much smaller than that in unsaturated soil. It can be assumed that there are two envelope curves in the equation of state of partially saturated soil if it is assumed that the volumetric residual deformation in partially saturated soil is mainly formed in the loading stage. One curve corresponds to the initial loading curve L_I, and the other curve corresponds to the complete compression curve L_C (see Fig. 1(b)). In the compression curve, the air content is smaller than that in the initial loading stage. A series of unloading–reloading curves L_U are included between the initial loading curve L_I and the fully compressed curve L_C. The expressions of the three curves are described in detail below.

• (1)

  Description of curve LI

Assuming that the volume compression relationship of the partially saturated soil skeleton part is

(3)

$$P_s=P_s\left({\mathrm{\Delta }}_s\right)=P_s\left(\mathrm{\Delta }\right)$$

The total volumetric compression deformation of each component in partially saturated soil is

(4)

$$\mathrm{\Delta }=\frac{\rho }{{\rho }_0}-1={\left({\alpha }_1\frac{{\rho }_{10}}{{\rho }_1}+{\alpha }_2\frac{{\rho }_{20}}{{\rho }_2}+{\alpha }_3\frac{{\rho }_{30}}{{\rho }_3}\right)}^{-1}-1$$

Where ρ₀ and ρ represent the density of soil before and after compression respectively. i = 1, 2, 3 represent air, water, and solid particles, respectively; α_i0 represents the percentage content of the i-th component; and ρ_i0 and ρ_i represent the density before and after compression, respectively, of the i-th component.

The relationship between the total hydrostatic pressure in the medium and the effective stress and pore hydrostatic pressure is

(5)

$$P=P_f+P_s$$

where the subscripts f and s represent the fluid and skeleton, respectively.

• (2)

  Description of curve L_C

According to Lyakhov, the volume of the skeleton in partially saturated soil is basically unchanged in the unloading stage under explosion loads. In this case,

(6)

$$P=P_f$$

• (3)

  Determination of curve L_U

Partially saturated soil is a kind of porous medium that can conveniently describe the load-unloading volume deformation relationship using the classic P-α model. In LS-DYNA (Zhao et al., 2004), the pore collapse model is a specific description of the P-α model. The model has been applied to the analysis of the ground shock load propagation law in TENSOR and has achieved good results.

The same assumption as the state equation described in LS-DYNA is adopted: the medium deformation is elastic when the volumetric strain follows Δ < Δ₁. When Δ > Δ₂, the pores in the partially saturated soil are completely closed, and the initial loading curve L_I and the complete compression curve L_C completely overlap. Physically, the bulk compressive modulus of the medium is much greater than the previous bulk compressive modulus, which means that the medium enters the approximate “locking” transformation stage. Between Δ₁ and Δ₂, the loading–unloading curve is determined by the following equation:

(7)

$$P_U(\mathrm{\Delta })=P_C\left(\frac{(1+{\mathrm{\Delta }}_B)(1+\mathrm{\Delta })}{1+{\mathrm{\Delta }}_{max}}-1\right)$$

where ${\mathrm{\Delta }}_B=P_C^{-1}(P_{max})$ represents the corresponding volume deformation when P_C = P_max in curve L_C, ${\mathrm{\Delta }}_{max}=P_I^{-1}(P_{max})$ represents the corresponding volume deformation when P_I = P_max in curve L_I, and the subscripts I, U, and C correspond to the curves L_I, L_U, and L_C, respectively.

The volume deformation can be expressed as

(8)

$$\overline{V}=\frac{V}{V_0}=\frac{{\rho }_0}{\rho }=\frac{1}{1+\mathrm{\Delta }}$$

Eq. (7) can be expressed as

(9)

$$P_U\left(\overline{V}\right)=P_C\left(\frac{{\overline{V}}_B}{{\overline{V}}_{min}}\overline{V}\right)$$ $${\overline{V}}_{min}=\frac{1}{{\mathrm{\Delta }}_{max}+1}$$

In the unloading–reloading curve L_U, the volume deformation modulus K_U can be expressed as

(10)

$$K_U=\frac{\partial P_U}{\partial \mathrm{\Delta }}\left(1+\mathrm{\Delta }\right)$$

For the unloading–reloading curve L_U,

(11)

$$\frac{\partial P_U}{\partial \mathrm{\Delta }}=\frac{\partial P_C({\mathrm{\Delta }}_a)}{\partial {\mathrm{\Delta }}_a}\frac{1+{\mathrm{\Delta }}_B}{1+{\mathrm{\Delta }}_{max}}$$

Therefore,

(12)

$$K_U=\frac{\partial P_C({\mathrm{\Delta }}_{\alpha })}{\partial {\mathrm{\Delta }}_a}\left(1+{\mathrm{\Delta }}_{\alpha }\right)$$ $${\mathrm{\Delta }}_a=\frac{(1+{\mathrm{\Delta }}_B)(1+\mathrm{\Delta })}{1+{\mathrm{\Delta }}_{max}}-1$$

In the actual calculation process, as long as the analytical description or the curve obtained from the data of the initial loading curve L_I and the fully compressed curve L_C is given, the unloading–reloading curve L_U can be directly calculated and determined by the program according to L_I and L_C.

(13)

$$P_I=P_f+P_S$$

For the P_f − Δ relationship,

(14)

$$\mathrm{\Delta }=1-\sum _{i=1}^3{\alpha }_i\frac{{\rho }_{i0}}{{\rho }_i}=1-\sum _{i=1}^3{\alpha }_i{\left(\frac{{\gamma }_i{P}_f}{{\rho }_{i0}{C}_{i0}^2}+1\right)}^{-\frac{1}{{\gamma }_i}}$$

(15)

$$\frac{d\mathrm{\Delta }}{d{P}_f}=\sum _{i=1}^3\frac{{\alpha }_i}{{\rho }_{i0}{C}_{i0}^2}{\left(\frac{{\gamma }_i{P}_f}{{\rho }_{i0}{C}_{i0}^2}+1\right)}^{-\frac{1}{{\gamma }_i}-1}$$

For the P_s − Δ relationship,

(16)

$$P_S=\int E_S(\mathrm{\Delta })d\mathrm{\Delta }$$

(17)

$$\frac{d{P}_S}{d\mathrm{\Delta }}=E_S$$

(18)

$$\frac{d{P}_I}{d\mathrm{\Delta }}=\frac{d{P}_f}{d\mathrm{\Delta }}+\frac{d{P}_S}{d\mathrm{\Delta }}$$

For curve L_C,

(19)

$$P_C=P_f$$

(20)

$$\mathrm{\Delta }=1-\sum _{i=1}^3{\alpha }_i\frac{{\rho }_{i0}}{{\rho }_i}=1-\sum _{i=1}^3{\alpha }_i{(\frac{{\gamma }_i{P}_C}{{\rho }_{i0}{C}_{i0}^2}+1)}^{-\frac{1}{{\gamma }_i}}$$

(21)

$$\frac{d\mathrm{\Delta }}{d{P}_C}=\sum _{i=1}^3\frac{{\alpha }_i}{{\rho }_{i0}{C}_{i0}^2}{(\frac{{\gamma }_i{P}_C}{{\rho }_{i0}{C}_{i0}^2}+1)}^{-\frac{1}{{\gamma }_i}-1}$$

The P-Δ relationship of L_U cannot be directly described by an analytical equation but can be solved numerically.

2.2

2.2 Validation of the P-α model

To validate the correctness of the model, LS-DYNA was used to carry out a numerical simulation of an experiment, which has been done before (Zhao et al., 2004). In the experiment, a partially saturated soil sample with a diameter of 0.9 m and a height of 4.74 m was made, and a Φ900 plane blast load generator is used to apply a uniform plane wave to the upper surface of the sample. Considered the axial symmetry of the sample in the experiment, the numerical model uses two-dimensional elements to reduce the size of elements and improve the calculation accuracy. The two-dimensional axisymmetric model with the side boundary as the normal constraint and the bottom boundary as the fixed boundary is established. The specific parameters of the numerical model and the P-α model are shown in Tables 2-4.

Table 2 Summary of the Physical Parameters of the Components.

	Density (kg/m3)	Wave velocity (m/s)	Isentropic index
Air	1.24	340	1.4
Fluid	1000	1500	3.0
Solid	2650	4500	7.0

Table 3 Summary of Component Contents at Various Depths.

Depth (m)	Air	Solid	Fluid
4.26 ∼ 4.74	0.020	0.570	0.410
3.79 ∼ 4.26	0.023	0.570	0.407
3.29 ∼ 3.79	0.026	0.570	0.404
2.8 ∼ 3.29	0.029	0.570	0.401
2.3 ∼ 2.8	0.032	0.570	0.398
1.8 ∼ 2.3	0.035	0.570	0.395
1.3 ∼ 1.8	0.038	0.570	0.392
0.8 ∼ 1.3	0.041	0.570	0.389
0.2 ∼ 0.8	0.044	0.570	0.386
0 ∼ 0.2	0.150	0.667	0.183

Table 4 Summary of the Mechanical Parameters of the Skeleton.

Parameter	υ	E1 (Pa)	E2 (Pa)	σy (Pa)
Value	0.3	1.0 × 107	0.6 × 107	1.0E5

The incident load in the numerical simulation is consistent with the pressure measured by the sensor on the soil sample surface in the experiment (Zhao et al., 2004). Fig. 2 and Fig. 3 show the incident load and validation of numerical simulation results, respectively. For the convenience of the following description, the vertical stress is simply referred to as the stress.

Close

Fig. 2

Stress-Time Curve of the Incident Load.

Export to PPT

Close

Fig. 3

Experimental and Numerical Simulation Results.

Export to PPT

For convenience, the first peak stress in Fig. 3 caused by only the incident load is called the incident peak stress. The second peak stress caused by the superposition of the incident load and the reflected load is called the reflected peak stress. The rise time at each depth in the numerical simulation is slightly less than the rise time obtained from the experiment at the same depth. The numerical data exhibit fluctuations in a certain range near the incident peak stress and the reflected peak stress. This phenomenon is caused by the bulk viscosity in the hydrodynamic equations, so the maximum stress of the fluctuation is not the actual stress at a certain depth (Von and Neumann r r., 1950). The intersection of the two straight lines between the rise stage and the gentle decrease stage after the peak stress is the actual peak stress at the depth. The numerical simulation results in Fig. 3 are not substantially different from the incident peak stress in the experimental results. The numerical simulation results are in good agreement with the experimental results for the decreasing trend and duration of the decreasing section after the incident peak stress. The above analysis shows that the waveform and peak stress of the numerical simulation results are in good agreement with the experimental results, which validates the reliability of the material model of the partially saturated soil material. Four numerical models with different air contents at the top were developed, in which the air content changed linearly with the depth in the model. Several loads with the same rise time, different peak stresses, and different duration times were applied to these models. Therefore, the propagation law of free-field plane shock waves in partially saturated soil can be studied.

3 Cases for numerical analysis

Assuming that the air content and other parameters of the soil are even in the horizontal direction, the simulatiton can be simplified to a one-dimensional problem under the action of a plane shock load such as a nuclear explosion. Based on the consideration of the material model of the partially saturated soil discussed before, the numerical model developed in this paper is a two-dimensional plane strain model. Considering the symmetry of the model, the size of the model in the horizontal direction is set to be much smaller than the size in the vertical direction to reduce the calculation time.

In the numerical simulation, two-dimensional plane strain models are established. The material model is the P-α model of the partially saturated soil, the solid content is set to 0.6, the sum of the air content and the water content is set to 0.4, the partially saturated soil is divided into many soil layers with equal thickness in the vertical direction, and the height of each thin layer is 0.01 m. The change in the air content of each layer is adjusted to simulate the phenomenon that the saturation of the partially saturated soil medium increases from top to bottom. In the numerical model, the air content of the partially saturated soil at the water-soil interface (depth y = 0 m) is denoted by α_y=0 and is selected according to the characteristics of the silty clay in the middle and lower reaches of the Yangtze River (Jiang and Luo, 2007), which are 0.01, 0.02, 0.03, and 0.04. The air content α decreases linearly in the depth direction, it is fully saturated at a depth of 15 m (y_α=0 = 15 m), and the bedrock is located at a depth of 20 m. Since the shock wave velocity in saturated soil can reach approximately 1600 m/s, the superposition of the reflected wave and the incident wave near the bedrock occurs in an extremely short period; therefore, the peak stress of the incident wave and the peak stress of the reflected wave cannot be accurately separated. A total of 81 elements are selected to record the stress time at various depths, and each element is separated by 0.2 m in the range of 0–18 m in depth. These stress-time curves are analysed so that the propagation law of shock waves under different geological conditions can be discussed.

The propagation law of shock waves in partially saturated soils with varying air contents is very different from that in homogeneous elastic media. Due to the increasing-hardening constitutive relation and the change in the propagation direction of partially saturated soil, the stress-time curve in partially saturated soil has a certain particularity. A total of three elements were selected in the model with α_y=0 = 0.01 under the shock wave with a peak stress of P_m = 0.2 MPa, rise time of t_r = 1 ms, and duration time of t_d = 200 ms. As shown in Fig. 4(a), the three elements are E1 at the top of the model; E2 at a depth of y = 6 m, which represents elements in the range of an air content greater than zero; and E3 at a depth of y = 18 m, which represents elements in the range of fully saturated soil. The stress-time curves of the three elements are shown in Fig. 4(b). The stress-time curve of E2 has fluctuations in a certain range near the incident peak stress. This phenomenon is caused by the bulk viscosity in the hydrodynamic equations, so the maximum stress of the fluctuation is not the actual stress at a certain depth. Since the stress of the incident load decreases linearly after reaching the peak stress, the stress at each depth tends to change linearly with time after the reaching the incident peak stress. Considering these reasons, a linear equation system, whose solution is taken as the incident peak stress at this depth, is established for the smooth decrease stage after the first fluctuation and the first rise stage. The reflected peak stress can be obtained by the same method. The incident peak stress and reflected peak stress are denoted by P₁ and P₂, respectively.

Close

Fig. 4

Typical Stress-Time Curves and Selected Elements.

Export to PPT

The stress-time curve of the top element E1 is almost the same as the stress-time curve of the incident load. The incident peak stress of element E2 is slightly larger than the peak stress of the incident load, reaching approximately 0.22 MPa, and the reflected peak stress reaches 0.41 MPa near 0.16 s, which exceeds the incident peak stress of the element. The incident peak stress of element E3 is approximately 0.73 MPa, which substantially exceeds the peak stress of the incident load. Since the element is only 2 m from the bedrock, the reflected load from the bedrock can be quickly attained. The reflected peak stress of element E3 reaches 1.4 MPa, and the stress decreases rapidly with time after the reflected peak stress. Based on a comparison of the incident peak stresses of the three elements, the incident peak stress gradually increases as the depth increases. The law of incident peak stress change with depth will be analysed in detail below.

4 Results and discussion

Incident loads with peak pressures of 0.2 MPa, 0.6 MPa, and 1.0 MPa are applied at the upper boundary of the model. The rise times of the incident loads are all 1 ms, and the duration times are 50 ms, 100 ms, 200 ms, and 300 ms. Four numerical models with different α_y=0 values under various incident loads are simulated by LS-DYNA. The incident peak stresses varying with depth in the four models are shown in Figs. 5-8.

Close

Fig. 5

Incident Peak Stress-Depth Curves in the Model with α_y=0 = 0.01.

Export to PPT

Close

Fig. 6

Incident Peak Stress-depth Curves in the Model with α_y=0 = 0.02.

Export to PPT

Close

Fig. 7

Incident Peak Stress-Depth Curves in the Model with α_y=0 = 0.03.

Export to PPT

Close

Fig. 8

Incident Peak Stress-Depth Curves in the Model with α_y=0 = 0.04.

Export to PPT

In Figs. 5-8, a similar relationship between the incident peak stress and depth in all numerical models with the same physical and geometric parameters under various incident loads is observed. In the range of 0–15 m, the incident peak stress decreases first and then increases with increasing depth or increases monotonically and reaches its maximum value at a depth of 15 m. This phenomenon is consistent with the results in (Zhao et al., 2004; Zhao et al., 2004). The incident peak stress barely changes as the depth increases at a depth greater than 15 m. The variation law of the incident peak stress in the range of 0–15 m is mainly controlled by two factors. The first factor is that the incident peak stress decreases with increasing propagation distance in partially saturated soil whose air content is constant (Von and Neumann r r., 1950), and the stress increases with depth caused by the increase in the wave impedance of the soil sample with depth (Lyakhov, 1974). When the decreasing rate of the incident peak stress with depth caused by the former is greater than the increasing speed of the stress caused by the latter, the incident peak stress gradually decreases with increasing depth. In contrast, the incident peak stress increases with increasing propagation distance.

In Figs. 5-8 (a)-(c), it is evident that the air content at the top of the partially saturated soil, α_y=0, has a certain influence on the variation law of the incident peak stress in the range of 0–15 m. Under the incident load with P_m = 1.0 MPa and t_d = 50 ms, the decrease rate and rise rate of the incident peak stress change in the range of 0–15 m both increase significantly as α_y=0 increases. The incident peak stress in the models with α_y=0 = 0.01 and α_y=0 = 0.02 under an incident load with P_m = 1.0 MPa and t_d = 100 ms in the range of 0–15 m increases exponentially as the depth increases. In models with α_y=0 = 0.03 and α_y=0 = 0.04, the incident peak stress decreases exponentially after a slight decrease in the range of 0–15 m. The law of the incident peak stress change discussed above is also applicable to models under incident loads with P_m = 0.2 MPa and P_m = 0.6 MPa. Regardless of the value of α_y=0, the incident peak stress in the range of 0–15 m increases exponentially as the depth increases in the models under incident loads with t_d = 200 ms and 300 ms. A larger α_y=0 leads to a greater increasing rate of the incident peak stress in the range of 0–15 m and a greater value of the maximum incident peak stress at 15 m. The comparison of the incident peak stress–depth curves from models with different α_y=0 values shows that under the same incident load, the larger incident peak stress occurs in the range of fully saturated soil in a model with a larger α_y=0.

The duration of the incident load has a certain influence on the variation law of the incident peak stress in the range of 0–15 m. In Figs. 5-8, for models under incident loads with t_d = 50 ms, the incident peak stress clearly decreases as the depth increases in the range of 0–7.5 m. The incident peak stress exhibits only a slight decrease in the range of shallow depths in the models under incident loads with t_d = 100 ms. The incident peak stress is obviously greater than that in the models under incident loads with t_d = 50 ms at each depth. The incident peak stress in the range of 0–15 m increases exponentially as the depth increases in the models under incident loads with t_d = 200 ms. The incident peak stress is greater than that in the models under incident loads with t_d = 100 ms at each depth. The variation law of the incident peak stress with depth in the models under incident loads with t_d = 300 ms is consistent with that in the models under incident loads with t_d = 200 ms. The incident peak stress is only slightly greater than that in the models under incident loads with t_d = 200 ms. In the range of depths greater than 15 m, where the is the soil is fully saturated, the incident peak stress is also affected by the duration time of the incident load: the incident peak stress increases as the duration time of the incident load increases, and the influence of the duration time on the incident peak stress decreases with the increase in duration. The above discussion shows that the value of the duration time of the incident load within a certain range (less than 200 ms) has a great influence on the incident peak stress at various depths, and the incident peak stress increases significantly at each depth as the duration time of the incident load increases. The influence of the duration of the incident loads on the incident peak stress is slight when the duration is greater than 200 ms.

The peak stress of the incident load P_m has a small influence on the variation law of the incident peak stress in the range of 0–15 m, but it directly affects the value of the peak stress at each depth. The ratio of the incident peak stress at a certain depth to the peak stress of the incident load should be considered.

To further study the influence of α_y=0, P_m, and t_d on the incident peak stress P₂, the dimensionless parameter overpressure coefficient at depth y₀ is defined as

where P_2y0 is the incident peak stress at depth y₀.

Because of the increase in the saturation and wave impedance of the soil sample with depth, the incident peak stress reaches the maximum in the fully saturated soil (Bloom, 2006; Wang, 2021). Obviously, it can be seen from Figs. 5-8 that the maximum overpressure coefficient is in the fully saturated soil, which is defined as

where P_2m represents the average value of the incident peak stress in fully saturated soil.

It is necessary to determine the relationship between the incident peak stress at shallow depths and the peak stress of the incident load for shallowly buried structures. The ratio of the incident peak stress at y = 0 to the peak stress of the incident load in Figs. 5-8 approaches 1. The variation in the incident peak stress with depth at shallow depths (0–7 m) is monotonous. The ratio of the incident peak stress at y = 7 m to the peak pressure of the incident load η_y=7 is plotted in Fig. 9, which can determine the influence of α_y=0, P_m, and t_d on the incident peak stress P₂.

Close

Fig. 9

η_y=7 in Various Models.

Export to PPT

Fig. 9 shows that η_y=7 decreases with the increase in α_y=0 under each incident load. According to the magnitude of the influence of the duration time of the incident load on the rate of decrease in η_y=7, the curves in Fig. 9 can be divided into four groups corresponding to t_d = 50 ms, 100 ms, 200 ms, and 300 ms. The slopes of the four groups of curves show that the rate of η_y=7 decreases with increasing t_d. This phenomenon indicates that the air content in partially saturated soil does not have a substantial effect on the incident peak stress in the range of shallow depths under long-duration loads such as nuclear explosions. However, when the incident load is a short-duration load, such as a chemical explosion, the air content in partially saturated soil has a great influence on the incident peak stress in the range of shallow depths.

The duration time of the incident load t_d can significantly affect the variation law of the incident peak stress in the range of shallow depths of the model under incident loads with a certain value of the peak stress. In the model with α_y=0 = 0.01, η_y=7 gradually increases as the incident load duration t_d increases, and the increasing rate decreases with increasing t_d. In the model with α_y=0 = 0.01 and under incident loads with t_d = 50 ms, η_y=7 is much smaller than that under incident loads with t_d = 100 ms. However, in the model with α_y=0 = 0.01 under incident loads with t_d = 300 ms, η_y=7 is not much greater than that under incident loads with t_d = 200 ms. The influence of t_d on η_y=7 discussed above is also appliable to models with different α_y=0 values, and the magnitude of the influence increases as α_y=0 increases. The above analysis shows that the duration time of the incident loads can significantly affect the incident peak stress variation law in the shallow depth range, and the influence gradually weakens as the duration increases.

Fig. 9 shows that, in the model under incident loads with a certain duration time, η_y=7 would increase slightly with increasing P_m, which is consistent with the existing experimental results.

Shock waves hardly attenuate in fully saturated soil (Schanz and Cheng, 2000); therefore, it is meaningful to study the ratio of the peak stress in fully saturated soil to the peak stress of the incident load for the underground structure in the fully saturated soil or below the fully saturated soil.

Fig. 10 shows that η_m increases with the increase in α_y=0 under each incident load. According to the magnitude of the influence of the duration time of the incident load on the rate of increase in η_m, the curves in Fig. 10 can be divided into four groups corresponding to t_d = 50 ms, 100 ms, 200 ms, and 300 ms. The slopes of the η_m curves for models under incident loads with t_d = 50 ms approach zero, indicating that the air content at the top of the models barely affects the incident peak stress in fully saturated soil. In the models under incident loads with t_d = 100 ms, η_m increases with increasing α_y=0, and the increasing rate gradually decreases. The trends of the curves of η_m in the models under incident loads with t_d = 200 ms and 300 ms do not vary substantially. Both trends exhibit a substantial increase with the increase in α_y=0, and there is a certain difference in the value of η_m. This shows that when the model is under a short-duration incident load such as a chemical explosion, the air content at the top of the soil does not substantially affect the incident peak stress in fully saturated soil at deep depths. However, the air content in partially saturated soil would significantly affect the incident peak stress in fully saturated soil in models under long-duration loads such as nuclear explosions. This law is exactly the opposite of the variation law of η_y=7 with α_y=0.

Close

Fig. 10

η_m in Various Models.

Export to PPT

The four groups of curves in Fig. 10 all show great influences of P_m, α_y=0 and t_d on η_m. η_m decreases significantly as P_m increases from 0.2 MPa to 0.6 MPa, and η_m exhibits a slight decrease as P_m increases from 0.6 MPa to 1.0 MPa. The value of α_y=0 of the models under incident loads with t_d = 50 ms has a slight influence on η_m. However, η_m increases substantially as α_y=0 increases in models under incident loads with t_d = 200 ms and 300 ms. In the model with α_y=0 = 0.01 under incident loads with t_d = 300 ms, η_m is not much greater than that in the same model under incident loads with t_d = 200 ms. However, the values of η_m in the model with α_y=0 = 0.04 increase significantly as the duration of the incident load changes from 300 ms to 200 ms. The increase in the peak stress of the incident load P_m causes η_m to decrease to a certain extent, but the incident peak stress in fully saturated soil generated by the incident load with a smaller P_m is still smaller than that generated by the incident load with a larger P_m, which is evident in Figs. 5-8.

The three variables, P_m, t_d, and α_y=0, all affect the maximum incident peak stress and the maximum overpressure coefficient in the fully saturated soil. The increase in α_y=0 will reduce the influence of P_m on η_m and increase the influence of t_d on η_m, and vice versa. Considering Figs. 5-8 and the effects of the above three factors on the incident peak stress in fully saturated soil, it is evident that underground structures will be at greater risk when the structure is located in partially saturated soil with a large value of α_y=0 or below fully saturated soil exposed to long-duration incident loads such as nuclear explosions.

Considering Figs. 9 and 10, long-duration incident loads are more threatening to underground structures than short-duration incident loads. Underground structures buried in fully saturated soil or below fully saturated soil with higher air contents above are more at risk than underground structures buried in partially saturated soil with smaller air contents.

5 Conclusions

In this paper, the gas content of partially saturated soil would decrease under a shock load was considered, which would cause volumetric plastic deformation, Lyakhov's constitutive model was improved. The P-α state equation was used to describe the volume compression relationship of partially saturated soils, and the numerical results were compared with experimental results to verify the correctness of the material model. Several partially saturated soil models in which the air content changes linearly in the depth direction were developed using the material model. The propagation law of incident loads with various peak stresses and durations in partially saturated soils was studied. The main conclusions are summarized as follows:

• (1)

  The peak stress and duration of the incident load and the air content of partially saturated soil all affect the ratio of the incident peak stress at various depths to the peak stress of the incident load.

• (2)

  The incident peak stress in partially saturated soil shows a trend of decreasing first and increasing later under an incident load with a shorter duration, and a trend of increasing monotonically under the influence of an incident load with a longer duration is observed. The rate of the incident peak stress change is proportional to the air content at the top of partially saturated soil. The incident peak stress reaches its maximum in fully saturated soil, which is 2.2 to 5.5 times the peak stress of the incident load.

• (3)

  Long-duration incident loads such as nuclear explosions are more threatening to underground structures than short-duration incident loads. Underground structures buried in fully saturated soil or below fully saturated soil with higher air contents above are more at risk than underground structures buried in partially saturated soil with smaller air contents.

In summary, the present study deals with the propagation law and influencing factors of plane shock waves in soil with varying saturation, which provides a certain basis for the related research on the dynamic response of underground structures in the future.