## 3 Numerical Solution Procedure
### The Governing Equation of the Soil Around the Pile Tip
Using the directional derivative rule (5), the governing equation along any curve coordinate \(S_{a+\kappa}S_{a-\kappa}\) at the computing node shown in Figure 6a can be written as Equation (6):  
\[\left\{ \begin{array}{l}\frac{\partial\xi}{\partial\tau}=\frac{1}{\sin 2\kappa}\Big{[}\sin(\alpha+\kappa)\frac{ \partial F}{\partial S_{a-\kappa}}-\sin(\alpha-\kappa)\frac{\partial F}{ \partial S_{a+\kappa}}\Big{]}\\ \frac{\partial F}{\partial z}=\frac{1}{\sin 2\kappa}\Big{[}-\cos(\alpha+ \kappa)\frac{\partial F}{\partial S_{a-\kappa}}+\cos(\alpha-\kappa)\frac{ \partial F}{\partial S_{a+\kappa}}\Big{]}\end{array} \right., \tag{5}\]  
\[\left\{ \begin{array}{l}\frac{\partial\sigma}{\partial S_{a+\kappa}}-\frac{1} {\sin 2\kappa}\left[\begin{array}{l}\frac{\partial R}{\partial S_{a-\kappa}}\sin(2 \psi-2\alpha-2\kappa)-\\ \frac{\partial R}{\partial S_{a+\kappa}}\sin(2\psi-2\alpha)\end{array}\right] -\frac{2R}{\sin 2\kappa}\left[\begin{array}{l}\frac{\partial\psi}{ \partial S_{a-\kappa}}\cos(2\psi-2\alpha-2\kappa)-\\ \frac{\partial\psi}{\partial S_{a+\kappa}}\cos(2\psi-2\alpha)\end{array}\right] \\ =f_{r}\cos(\alpha+\kappa)+f_{z}\sin(\alpha+\kappa)\\ \frac{\partial\sigma}{\partial S_{a-\kappa}}-\frac{1}{\sin 2\kappa}\left[\begin{array}{l} \frac{\partial R}{\partial S_{a-\kappa}}\sin(2\psi-2\alpha)-\\ \frac{\partial R}{\partial S_{a+\kappa}}\sin(2\psi-2\alpha+2\kappa)\end{array} \right]-\frac{2R}{\sin 2\kappa}\left[\begin{array}{l}\frac{\partial\psi}{ \partial S_{a-\kappa}}\cos(2\psi-2\alpha)-\\ \frac{\partial\psi}{\partial S_{a+\kappa}}\cos(2\psi-2\alpha+2\kappa)\end{array} \right]\\ =f_{r}\cos(\alpha-\kappa)+f_{z}\sin(\alpha-\kappa)\end{array} \right., \tag{6}\]  
where \(\alpha\) is the angle between the angular bisector and the positive direction of the \(r\)-axis, \(\kappa\) is the angle between the curve coordinate direction and the angular bisector, \(f_{r}=-(R/r)\cdot(1+\cos 2\psi)\), and \(f_{z}=\gamma-(R/r)\cdot\sin 2\psi\). The governing Equation (6) describes the variation of stress state variables \((\sigma,R,\psi)\) along curve coordinate direction \(S_{a+\kappa}\) and \(S_{a-\kappa}\). To calculate the position and stress information \((\sigma,R,\psi,r,z)\) of an unknown node from the known nodes, the discrete form of the governing equation is then derived. For the known nodes A and B and node C to be calculated as shown in Figure 6b, the midpoint  
Figure 6: (**a**) Curve coordinate system at computing node A; (**b**) discrete computing nodes. M1 of AC in the curve coordinate \(S_{A_{A}+\kappa_{A}}S_{A_{A}-\kappa_{A}}\) and the midpoint M2 of BC in the curve coordinate \(S_{A_{B}+\kappa_{B}}S_{A_{B}-\kappa_{B}}\) are selected to make a center difference scheme. The difference formulas are shown in Equations (7) and (8):  
\[\frac{\partial\sigma}{\partial S_{\alpha_{M_{1}}+\kappa_{M_{1}}}}=\frac{ \sigma_{\mathrm{C}}-\sigma_{A}}{\Delta AC}\,\ \frac{\partial R}{\partial S_{\alpha_{M_{1}}+\kappa_{M_{1}}}}=\frac{R_{ \mathrm{C}}-R_{A}}{\Delta AC}\,\ \frac{\partial\psi}{\partial S_{\alpha_{M_{1}}+\kappa_{M_{1}}}}=\frac{ \psi_{\mathrm{C}}-\psi_{A}}{\Delta AC}, \tag{7}\]  
\[\frac{\partial\sigma}{\partial S_{\alpha_{M_{2}}-\kappa_{M_{2}}}}=\frac{ \sigma_{\mathrm{C}}-\sigma_{B}}{\Delta BC}\,\ \frac{\partial R}{\partial S_{\alpha_{M_{2}}-\kappa_{M_{2}}}}=\frac{R_{ \mathrm{C}}-R_{B}}{\Delta BC}\,\ \frac{\partial\psi}{\partial S_{\alpha_{M_{2}}-\kappa_{M_{2}}}}=\frac{ \psi_{\mathrm{C}}-\psi_{B}}{\Delta BC}, \tag{8}\]  
Substituting Equations (7) and (8) in (6), the discrete governing equation is obtained as follows:  
\[\left\{\begin{array}{l}\frac{\sigma_{\mathrm{C}}-\sigma_{A}}{\Delta AC}- \frac{1}{\sin 2\kappa_{M_{1}}}\cdot\left[ \begin{array}{c}\frac{R_{ \mathrm{C}}-R_{B}}{\Delta BC}\sin\bigl{(}2\psi_{M_{1}}-2\alpha_{M_{1}}-2\kappa _{M_{1}}\bigr{)}-\\ \frac{R_{\mathrm{C}}-R_{A}}{\Delta AC}\sin\bigl{(}2\psi_{M_{1}}-2\alpha_{M_{1}} \bigr{)}\end{array} \right]+\frac{2R_{M_{1}}}{\sin 2\kappa_{M_{1}}}.\\ \\ \left[