## 3 Numerical Solution Procedure  
### The Governing Equation of the Soil Around the Pile Tip  
The failure zone below the pile end plane is determined by the zero-extension line combined with the assumed failure pattern mentioned above and the boundary conditions. Investigating the end-bearing capacity of a single pile shown in Figure 3, the governing equation of soil element around the pile tip in the cylindrical coordinate system \((r,\theta,z)\) can be written as follows:  
\[\left\{ \begin{array}{l}\frac{\partial\tau_{r}}{\partial r}+\frac{\partial \tau_{r}}{\partial z}+\frac{\sigma_{r}-\sigma_{\theta}}{r}=0\\ \frac{\partial\tau_{r}}{\partial r}+\frac{\partial\tau_{r}}{\partial z}+\frac{ \tau_{r}}{r}-\gamma=0\end{array} \right., \tag{1}\]  
where \(\sigma_{r}\), \(\sigma_{z}\), \(\sigma_{\theta}\), and \(\tau_{rz}\) are components of the stresses; and \(\gamma\) is the soil unit weight. Figure 3: The cylindrical coordinate system and the analyzed soil element. Figure 2: (**a**) The division of soil around the pile tip; (**b**) the load acting on the failure zone boundary. The state of stress at a point is shown in Figures 4 and 5. Using the Mohr-Coulomb yield criterion, the components of the stresses can be written as follows, with the assumption of \(\sigma_{\theta}\) being equal to the minor principal stress [24]:  
\[\left\{ \begin{array}{l}\sigma_{z}=\sigma-R\cos 2\psi=\sigma(1-\sin\phi\cos 2 \psi)-c\cos\phi\cos 2\psi\\ \sigma_{r}=\sigma+R\cos 2\psi=\sigma(1+\sin\phi\cos 2\psi)+c\cos\phi\cos 2 \psi\\ \tau_{rz}=R\sin 2\psi=\tau_{zr}=\sigma\sin\phi\sin 2\psi+c\cos\phi\sin 2 \psi\\ \sigma_{\theta}=\sigma_{3}=\sigma-R=\sigma(1-\sin\phi)-c\cos\phi\end{array}  \right., \tag{2}\]  
where \(\sigma\) is the mean stress (\(\sigma=(\sigma_{r}+\sigma_{z})/2=(\sigma_{1}+\sigma_{3})/2\)), \(R\) is the radius of the Mohr stress circle (\(R=(\sigma_{1}-\sigma_{3})/2\)), \(\psi\) is the angle between the direction of the major principal stress and the positive direction of the \(r\)-axis, \(\phi\) is the mobilized friction angle, and \(c\) is the soil cohesion intercept. Substituting Equation (2) in Equation (1), the governing equation expressed in terms of \(\sigma\) and \(R\) is obtained as follows:  
\[\left\{ \begin{array}{l}\frac{\partial\sigma}{\partial\tau}+\frac{\partial R }{\partial\tau}\cos 2\psi+\frac{\partial R}{\partial\varepsilon}\sin 2\psi+2R \Bigl{[}-\sin 2\psi\frac{\partial\psi}{\partial\tau}+\cos 2\psi\frac{ \partial\psi}{\partial\varepsilon}\Bigr{]}=-\frac{R}{r}(1+\cos 2\psi)\\ \frac{\partial\sigma}{\partial\varepsilon}+\frac{\partial R}{\partial\tau} \sin 2\psi-\frac{\partial R}{\partial\varepsilon}\cos 2\psi+2R\Bigl{[}\cos 2\psi \frac{\partial\psi}{\partial\tau}+\sin 2\psi\frac{\partial\psi}{\partial\varepsilon} \Bigr{]}=\gamma-\frac{R}{r}\sin 2\psi\end{array} \right., \tag{3}\]  
Figure 4: The state of stress at a point and the directions of the stress characteristic line and zero-extension line (ZEL). Figure 5: Mohr stress circle at intermediate state and ultimate state. The zero-extension line along which the linear incremental strain is zero is defined by two families that intersect with the direction of major principal stress at \(\xi\), as shown in Figure 4. The equations for the zero-extension line can be written as follows:  
\[\left\{ \begin{array}{l}+ZEL\ direction:\frac{dz}{dr}=\tan(\psi+\zeta)\\ -ZEL\ direction:\frac{dz}{dr}=\tan(\psi-\zeta)\end{array} \right., \tag{4}\]  
where \(\xi=\pi/4-\
u/2\) and \(v\) is the dilation angle. The ZEL net can be constructed by simultaneously solving the soil governing Equation (3) and the ZEL Equation (4) from the boundary between zone I and zone II to the boundary between the rigid cone and zone I. The finite-difference method is used to solve these equations and the difference forms of these equations are derived to obtain the position and stress information of the node in the failure zone.