Ultimate bearing capacity of circular shallow foundations in frozen clay

3.1. Test device and procedure [29]

The model test system mainly includes refrigeration system, circulation system, data acquisition system, control system and test platform etc. (as shown in Fig. 1). The experimental platform is a circular pit in the Institute of Frozen Soils of Anhui University of Science and Technology. Its radius and depth are 1.5 m and 0.9 m, respectively. The soil samples are Huainan clay, and its basic physical parameters are listed in Table 1. A high-precision refrigeration system using alcohol was used in the experiment, and the soil was reduced to the specified temperature by artificial freezing tubes. The temperature of the freezing tubes was –15 °C, and the layout of the freezing tubes are shown in Fig. 2. The circular rigid plate whose load was directly applied to the frozen soil by the hydraulic jack, and the fast load-keeping method was adopted, which was divided into 8 steps with an interval of 1 h. The data acquisition used a high-precision YE2532 static strain gauge, and the whole process of the test was controlled by computer. For comparison, we also carried out a model test of the foundation bearing capacity of clay, and the specific arrangement is consistent with the frozen soil.

Fig. 1Model test system

Fig. 2Layout of freezing tubes

3.2. Test results

The load-settlement (P-S) curve of frozen clay at –15°C and conventional clay are shown in Fig. 3 and 4. As we can see:

(1) The strength of the frozen clay is larger and the ultimate bearing capacity is stronger than unfrozen soil, which is due to the water in the soil are converted into ice at low temperature. At this time, the strength is mainly affected by the cementation of the ice, which is also the main influence of temperature on the strength of frozen soil.

(2) Those load-settlement curves of frozen clay all have proportional limitation of 170 kN and are equal in size, this shows that the ultimate bearing capacity of frozen clay is about 2.41 MPa.

(3) Conventional clay fails by general shear mechanisms, and frozen clay fails by local shear mechanisms. This could be determined by the P-S curves, which the P-S curve of local shear mechanisms often haven’t obvious turning point, and it’s mainly characterized by deformation, the sliding surface did not develop to the ground.

Fig. 3P-S curve of frozen clay at –15°C

Fig. 4P-S curve of conventional clay

4.1. The generalized Kelvin model

The creep characteristic of frozen soil is complex. Many scholars have studied it, and several frozen soil creep models have been established [30-34]. The creep curves for frozen clay can be described by the generalized Kelvin model consisting of a Hooke element and a Kelvin elements (Fig. 5).

Fig. 5The generalized Kelvin model

The constitutive equation is:

The creep equation is:

where ${\epsilon }^e$ and ${\epsilon }^{ve}$ represent the elastic strain and the viscoelastic strain, respectively; $E_1$ is the elasticity parameter of spring in the generalized Kelvin body; $E_k$ is the elasticity parameter of spring in Kelvin body; $\eta$ is the viscosity coefficient of a dashpot in the generalized Kelvin body; $t$ is time.

4.2. Model flexibility matrix [35]

It can be seen from Fig. 5 that strain of the generalized Kelvin model can be decomposed into elastic strain and viscoelastic strain. Elastic strain can be determined by the generalized Hooke’s law. For three-dimensional axisymmetric problems, assuming the stress boundary conditions remain unchanging, and the Poisson’s ratio does not change with time [36], the elastic strain can be expressed as:

Assuming that the soil is isotropic, the three-dimensional the matrix [$C$] can be express:

Considering the rheological properties of the generalized Kelvin model, the basic form of the creep equation is obtained from the rheological constitutive equation by substituting the differential operator and solving the differential equation:

Rewrite the Eq. (17) and (19) into tensor form:

Thus, the three-dimensional differential constitutive relationship of the generalized Kelvin model under the constant Poisson's ratio can be expressed as:

4.3. Creep test of frozen clay

The uniaxial creep tests of frozen clay were graded loading, and the loading coefficients were 0.3, 0.5, and 0.7, respectively. The specimen was a cylinder with a diameter of 50 mm and a height of 100 mm. It was prepared by multi-layer wet tamping method (as shown in Fig. 6). After demoulded, it was placed in a low temperature box at –15°C for 48 hours and taken out for testing. Fig. 7 shows a damage specimen after the test. The test results are shown in Fig. 8.

Fig. 6Specimen in creep test

Fig. 7Specimen after test

Simulated annealing algorithm was used to optimize the parameters of the general Kelvin model [37]. The specific parameters are shown in Table 2. The uniaxial creep test results and calculated results of the generalized Kelvin model of artificially frozen clay at –15 °C are shown in Fig. 8. It can be seen that the generalized Kelvin model could well describe the creep of frozen clay.

Fig. 8Experimental and calculated results of uniaxial creep tests of artificially frozen clay at –15 °C

5.1. Numerical model

The finite element program ANASY is applied for the numerical simulation. The finite element model is consistent with the model test, that are its radius and depth are 1.5 m and 0.9 m, respectively. The eight-nodes hexahedral element type with independent temperature degrees of freedom are adopted and it has plasticity, creep, swelling, stress stiffening, large deflection and large strain capabilities. The mesh size is between 140 mm×70 mm and 20 mm×20 mm. After meshing, there are 10669 elements generated. The FEM model and mesh division are shown in Fig. 9.

Fig. 9Finite element meshes

a)

b)

5.2. Parameters

The parameters of the soil, such as thermal conductivity and specific heat, are changing with the state of the soil. Consequently, the influence of the phase transformation on the soil parameters should be taken into considered. The parameters are shown in Table 3.

5.3. Numerical test results

Firstly, a three-dimensional solid model was established, and the temperature load was applied to the model. Then the step loading method was applied to simulate the multi-stage loading, and loading mode was consistent with the model test loading.

The comparison between simulation results and model test results are shown in Fig. 10. As can be seen that the simulation results are consistent with model test results. Those curves are composed of three stages, but the settlement values have a certain gap. When the load is small, numerical results are smaller than experimental results, this is mainly due to the foundation soil produces initial settlement under the action of its own weight, but the numerical model don’t take initial settlement into account; as load increases, numerical results gradually greater than experimental results, the frozen soil is regarded as ideal elastoplastic body when simulates. But in fact, the properties of the frozen soil are very complicated, which is closer to viscoelastic material and difficult to fully describe its properties with a simple component model.

Meanwhile, the generalized Kelvin model often overestimates the third-stage creep of the soil. In this paper, we pay more attention to the failure mechanism and ultimate bearing capacity of frozen soil foundation, the creep properties of the frozen soil can be approximated describe by the generalized Kelvin model. And the room temperature is higher than the freezing temperature during testing. Although we had taken some measures to maintain the heat, there is inevitably heat transfer between the soil and air, which brings difficulties to our numerical simulation, so it is simplified, this also leads to some differences in the results of the simulation. Compared with the P-S curves, there are obvious proportional boundary points. The boundary is bounded by the point, and both segments are close to linear, which is consistent with the P-S curve law of local shear failure, this is consistent with the observation from model test. In addition, the compaction and freezing degree of the soil in model test can also affect the result. The bearing capacity can be determined by the proportional limit of the P-S curves, the numerical result is about 190 KN, which is larger than experimental result of 170 KN. The error determined by the two methods is within 15 %, therefore, it could be considered that the simulation results are in agreement with the actual.

Fig. 10Comparison between simulation values and model test values

6.1. Basic assumptions [17, 21, 24, 38]

(1) The soil is weightless media.

(2) The base in foundations is completely smooth.

(3) The failure plane of local shear failure conforms to the Prandtl sliding failure surface, and the sliding zone consists of the active zone I, the transition zone II and the passive zone III. It is shown in Fig. 11, and $d$ is the radius of circular foundation:

where $\phi$ is the internal friction angle in soil.

(4) In the range of embedding depth of foundation, the friction between the foundation and soil and the shear strength of the soil on both sides of the foundation are not taken into account, and the weight of the buried soil on both sides is simplified as overloading.

Fig. 11Idealized failure mechanism.

6.2. Analytical solution

6.2.1. Analysis of active zone

For the circular foundation, under the assumption of the local shear failure mode, the force of the active zone is shown in Fig. 12. The positive stress on the active failure surface is ${\sigma }_a$, the shear stress is ${\tau }_a$, and it is equal to the cohesion.

Fig. 12External forces acting on the active area

As shown in the Fig. 12, the surface area of the cone (except the area of the top) is:

23

$S_1=\pi dd\mathrm{s}\mathrm{e}\mathrm{c}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)=\pi d^2\mathrm{s}\mathrm{e}\mathrm{c}\left(\frac{\pi }{4}+\frac{\phi }{2}\right).$

From the vertical force balance condition of the active zone, it is obtained:

24

$P_u\pi d^2={\sigma }_a{S}_0\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)+{\tau }_a{S}_0\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)=S_0\left[{\sigma }_a\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)+C\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)\right],$

that is:

25

$P_u=\mathrm{s}\mathrm{e}\mathrm{c}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)\left[{\sigma }_a\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)+C\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)\right]={\sigma }_a+C\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right).$

6.2.2. Analysis of transition zone and passive zone

Since the local shear failure fracture surface extends into the passive zone, it can be considered conservatively that the fracture plane only develops to the boundary between the transition zone and the passive zone, that is, point B in Fig. 11. Compared with the general shear failure mode, the main difference between the two modes is that the pressure in the BDD₁B₁ plane changes from passive earth pressure $P_p$ to non-limiting Rankine passive earth pressure $P_p^{*}$ [39]. A micro-body of the transition zone and the passive zone is subjected to force analysis (Fig. 13). In the figure, BC is a log spirals, where $r$ is the radius of the spiral, $r=AC\mathrm{e}\mathrm{x}\mathrm{p}\left(\varphi \mathrm{t}\mathrm{a}\mathrm{n}\phi \right)$. The moment balance conditions of the micro-body to point A are:

(1) The area of plain ADD₁A₁ is:

26

$S_1=\frac{\left(A{A}_1+D{D}_1\right)}{2}AD=\frac{\left[2+\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\frac{\pi }{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi }\right]}{2}{d}^{2}d\theta \mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\frac{\pi }{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi }.$

The force on plain ADD₁A₁ is the ground overload $q$, and the moment on point A is:

27

$\begin{array}{c}{M}_1=q{S}_1\frac{A{A}_1+2D{D}_1}{3\left(A{A}_1+D{D}_1\right)}AD=\frac{1}{6}q{d}^3\mathrm{t}\mathrm{a}{\mathrm{n}}^2\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\pi \mathrm{t}\mathrm{a}\mathrm{n}\phi }\left[3+2\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\frac{\pi }{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi }\right]d\theta .\end{array}$

(2) For the area of plain BD₁DB₁, the length of BB₁ can be approximated to be equal to the length of DD₁:

28

$S_2=D{D}_{1}BD=d^2\left(1+\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\frac{\pi }{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi }\right)d\theta e^{\frac{\mathrm{\pi }}{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi }.$

The force on the plain BD₁DB₁ is the non-limiting Rankine passive earth pressure, and its specific expression is:

29

$\begin{array}{c}{P}_p^{\mathrm{*}}=q{K}_0+\frac{2z\left[q\left(K_p-K_0\right)+2c\sqrt{K_p}\right]}{\left(d{e}^{\frac{\pi }{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi }\right)},K_p=\mathrm{t}\mathrm{a}{\mathrm{n}}^2\left(\frac{\pi }{4}+\frac{\phi }{2}\right),K_0=1-\mathrm{s}\mathrm{i}\mathrm{n}\phi .\end{array}$

The moment of force on A is:

30

$\begin{array}{l}{M}_2={\int }_0^{BD}\left[\frac{l}{BD}\left(q{K}_p+2c\sqrt{K_p}-q{K}_0\right)+q{K}_0\right]D{D}_{1}dl\cdot l=\frac{1}{6}{d}^3{e}^{\pi \mathrm{t}\mathrm{a}\mathrm{n}\phi }d\theta \\ ·\left[2q\mathrm{t}\mathrm{a}{\mathrm{n}}^2\left(\frac{\pi }{4}+\frac{\phi }{2}\right)+4c\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)+q\left(1-\mathrm{s}\mathrm{i}\mathrm{n}\phi \right)\right]\left[1+\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\frac{\pi }{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi }\right].\end{array}$

(3) The force on the surface CBB₁ is normal stress ${\sigma }_t$, the direction of action points to the pole A of the helix, and the shear stress ${\tau }_t$, the magnitude is equal to the cohesion, and is evenly distributed along the boundary. The moment of the surface force on the point A is:

31

$\begin{array}{l}{M}_3={\int }_0^{\varphi }d\mathrm{s}\mathrm{e}\mathrm{c}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\beta \mathrm{t}\mathrm{a}\mathrm{n}\phi }d\beta d\theta {\int }_0^{\frac{\pi }{2}}{d}^2\mathrm{s}\mathrm{e}{\mathrm{c}}^2\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{2\varphi \mathrm{t}\mathrm{a}\mathrm{n}\phi }\\ \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}·d\varphi =c{d}^3\mathrm{s}\mathrm{e}{\mathrm{c}}^3\left(\frac{\pi }{4}+\frac{\phi }{2}\right)d\theta \mathrm{c}\mathrm{o}{\mathrm{t}}^3\phi \left[\frac{1}{3}{e}^{\frac{3}{2}\pi \mathrm{t}\mathrm{a}\mathrm{n}\phi }-\frac{1}{2}{e}^{\pi \mathrm{t}\mathrm{a}\mathrm{n}\phi }+\frac{1}{6}\right].\end{array}$

(4) The area of plain CAA₁ is:

32

$S_4=\frac{1}{2}A{C}^{2}d\theta =\frac{1}{2}{d}^2\mathrm{s}\mathrm{e}{\mathrm{c}}^2\left(\frac{\pi }{4}+\frac{\phi }{2}\right)d\theta .$

The moment of force on A is:

33

$M_4={\sigma }_a{S}_4\frac{1}{3}AC=\frac{1}{6}{d}^3\mathrm{s}\mathrm{e}{\mathrm{c}}^3\left(\frac{\pi }{4}+\frac{\phi }{2}\right)d\theta {\sigma }_a.$

Due to static equilibrium conditions $\sum M_A=\text{0}$, there are:

34

$M_4=M_1+M_2+M_3.$

Substituting Eqs. (27), (30), (31) and (33) into Eq. (34), as follows:

35

$\begin{array}{l}{\sigma }_a=c\mathrm{c}\mathrm{o}{\mathrm{t}}^3\phi \left(2e^{\frac{3}{2}\pi \mathrm{t}\mathrm{a}\mathrm{n}\phi }-3e^{\pi \mathrm{t}\mathrm{a}\mathrm{n}\phi }+1\right)+c4\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)\mathrm{c}\mathrm{o}{\mathrm{s}}^2\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\pi \mathrm{t}\mathrm{a}\mathrm{n}\phi }\\ ·\left[1+\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\frac{\pi }{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi }\right]+q\mathrm{s}\mathrm{i}{\mathrm{n}}^2\left(\frac{\pi }{4}+\frac{\phi }{2}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)\left[3+2\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\frac{\pi }{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi }\right]\\ ·e^{\pi \mathrm{t}\mathrm{a}\mathrm{n}\phi }+q\left[2\mathrm{t}\mathrm{a}{\mathrm{n}}^2\left(\frac{\pi }{4}+\frac{\phi }{2}\right)+1-\mathrm{s}\mathrm{i}\mathrm{n}\phi \right]\mathrm{c}\mathrm{o}{\mathrm{s}}^3\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\pi tan\phi }\\ ·\left[1+\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\frac{\pi }{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi }\right].\end{array}$

Substituting Eq. (35) into Eq. (25), for brevity, it can be written as:

36

$P_k=q{N}_q+c{N}_c.$

In which:

37

$N_q=\mathrm{s}\mathrm{i}{\mathrm{n}}^2\left(\frac{\pi }{4}+\frac{\phi }{2}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\pi \mathrm{t}\mathrm{a}\mathrm{n}\phi }\left[3+2\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\frac{\pi }{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi }\right]+e^{\pi \mathrm{t}\mathrm{a}\mathrm{n}\phi }\mathrm{c}\mathrm{o}{\mathrm{s}}^3\left(\frac{\pi }{4}+\frac{\phi }{2}\right)$
$\begin{array}{l}\\ \cdot \left[1+\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\frac{\pi }{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi }\right]\left[2\mathrm{t}\mathrm{a}{\mathrm{n}}^2\left(\frac{\pi }{4}+\frac{\phi }{2}\right)+1-\mathrm{s}\mathrm{i}\mathrm{n}\phi \right],\end{array}$

38

$N_c=\mathrm{c}\mathrm{o}{\mathrm{t}}^3\phi \left(2e^{\frac{3}{2}\pi \mathrm{t}\mathrm{a}\mathrm{n}\phi }-3e^{\pi \mathrm{t}\mathrm{a}\mathrm{n}\phi }+1\right)+4e^{\pi \mathrm{t}\mathrm{a}\mathrm{n}\phi }\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)\mathrm{c}\mathrm{o}{\mathrm{s}}^2\left(\frac{\pi }{4}+\frac{\phi }{2}\right)$
$·\left[1+\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right)e^{\frac{\pi }{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi }\right]+\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\pi }{4}+\frac{\phi }{2}\right).$

Fig. 13External forces acting on the transition zone and passive zone

6.2.3. Comparison with empirical formulas

Fig. 14 compares the results of $N_c$ and $N_q$ from Terzaghi empirical formula, Vesic empirical formula, and our method in this paper. As can be seen:

(1) The foundation bearing capacity coefficients $N_c$ and $N_q$ obtained by the proposed method increase with the increase of the internal friction angle, and the variation law are consistent with the Terzaghi and Vesic empirical formulas.

(2) The coefficients of circular foundation obtained by Terzaghi empirical formula is too small, which seems to be too conservative.

(3) The $N_c$ and $N_q$ obtained by the Vesic empirical method are closely related to $I_r$/$I_{r\left(cr\right)}$. The $N_q$ calculated by this paper is between 0.5 ≤ $I_r$/$I_{r\left(cr\right)}$ < 0.9, and $N_c$ is between 0.3 ≤ $I_r$/$I_{r\left(cr\right)}$ < 0.5.

Fig. 14Comparisons between Nq and Nc by different methods of circular foundation subjected to local shear failure mechanisms

a)$N_q$

b)$N_c$

In addition, Terzaghi’s and Vesic’s analytical solutions belong to empirical formulas, and the method proposed in this paper is derived from rigorous deduction under certain assumptions and has certain theoretical basis. Further studies are needed to take into account the gravity and friction of foundation.

Assume that the foundation is a weightless medium, the results of different methods are listed in Table 4. It can be seen that the bearing capacity of the circular foundation under the local shear failure mode obtained by different methods is quite different. The calculation results of the equation in this paper are between the model test results and the numerical simulation results, which verifies the rationality of the method in this paper.

Table 4Comparison of ultimate bearing capacity by different methods

Methods	Bearing capacity (MPa)
The expanding cavity model	2.906
Terzaghi empirical formula	1.659
Vesic empirical formula	3.223
The test in this paper	2.406
The FEM in this paper	2.689
The formula in this paper	2.568