## 2 Model Formulation and Analysis
### Mathematical Modeling
The vertical displacement of the surrounding soil \(W\) is expressed by the vertical surface displacement \(w(r)\) and the shape function \(\phi(z)\) as \[W(r,z)=w(r)\phi(z)\] (6) where the boundary conditions can be assumed as \(\phi(0)=1\) and \(\phi(H)=0\). Inserting Equations (6) and (4) into Equation (3) yields the following equation:  
\[ \begin{split}&\prod=\Pi_{c}+\int_{\tau_{c}}^{\frac{d}{2}}\int_ {0}^{2\pi}\int_{0}^{H}\left[\begin{array}{c}\frac{(1-\
u_{c})E_{c}}{2(1+ \
u_{c})(1-2\
u_{c})}\left(\frac{d\phi}{dx}\right)^{2}w^{2}\\ +\frac{E_{c}}{4(1+\
u_{c})}\phi^{2}\left(\frac{d\psi}{dr}\right)^{2}\end{array} \right]dzd\theta rdr\\ &-\int_{0}^{r_{c}}\int_{0}^{2\pi}\sigma_{0}wd\theta rdr\end{split}  \tag{7}\]  
By minimizing the function \(\Pi\) with respect to \(w\) based on the variational principles, the following governing equation is obtained. \[kw-G(\frac{d^{2}w}{dr^{2}}+\frac{1}{r}\frac{dw}{dr})=\sigma_{o}\quad r_{c} \leq r\leq\frac{d}{2} \tag{8}\]  
where  
\[ \begin{cases}k=\frac{(1-\
u_{c})E_{c}}{(1+\
u_{c})(1-2\
u_{c})}\int\limits_{0} ^{H}\left(\frac{d\phi}{dx}\right)^{2}dz\\ \quad\quad G=\frac{E_{c}}{2(1+\
u_{c})}\int\limits_{0}^{H}\phi^{2}dz\end{cases}  \tag{9}\]  
By minimizing the total potential energy \(\Pi\) of Equation (7) with respect to \(\phi\), the following equation is obtained:  
\[-m\frac{d^{2}\phi}{dz^{2}}+n\phi=0 \tag{10}\]where  
\[\left\{ \begin{array}{l}m=\frac{(1-u_{\rm s})E_{\rm s}}{(1+v_{\rm s})(1-2v_{\rm s} )}\int_{r_{\rm c}}^{\frac{\phi}{2}}w^{2}rdr\\ n=\frac{E_{\rm s}}{2(1+v_{\rm s})}\int_{r_{\rm c}}^{\frac{\phi}{2}}\Big{(}\frac {dw}{dr}\Big{)}^{2}rdr\end{array} \right. \tag{11}\]  
By introducing a new parameter \(\eta\), the following equation can be obtained:  
\[\frac{n}{m}=\Big{(}\frac{\eta}{H}\Big{)}^{2} \tag{12}\]  
Substitute Equation (12) into Equation (10), and the following equation can be derived:  
\[\frac{d^{2}\phi}{dz^{2}}-\Big{(}\frac{\eta}{H}\Big{)}^{2}\phi=0 \tag{13}\]  
Equation (11) implies that \(m\) and \(n\) depend on the surface displacement function \(\mathbf{w}\). At the same time, \(w\) depends on the variable \(\eta\). Therefore, the surface function \(\mathbf{w}\) and shape function \(\mathbf{\phi}\) are coupling together. To uncouple these functions, [PERSON] et al. [38] proposed an iterative method to solve Equations (8) and (12) simultaneously. [PERSON] et al. [39] solved Equation (13) with the boundary conditions \(\phi(0)=1\) and \(\phi(H)=0\), then revealed that the solution is the shape function in Equation (14). \[\phi(z)=\frac{\sinh[\eta(1-z/H)]}{\sinh(\eta)} \tag{14}\]  
Equation (14) is a general form for \(\phi(z)\), and [PERSON] et al. [40] proposed an approximate form of the shape function \(\phi(\rm z)\) in a thin layer as follows:  
\[\phi(z)=1-z/H \tag{15}\]  
According to Equation (8), the function, \(w(r)\), is the solution for a second-order variable coefficient differential equation, which consists of two parts, namely the homogeneous solution and the specific solution. The homogeneous solution is given by the Bessel function, which can be readily obtained. For the specific solution, the constant \(\frac{\phi_{\rm c}}{K}\) satisfies this equation. Therefore, its general solution is derived as  
\[\left\{ \begin{array}{l}w(r)=C_{1}K_{0}(a\times r)+C_{2}I_{0}(a\times r)+ \frac{\sigma_{\rm c}}{K}\\ a=\sqrt{\frac{k}{C}}\end{array} \right. \tag{16}\]  
where \(I_{0}\) and \(K_{0}\) are the zero-order modified Bessel functions of the first and the second kind, respectively. Only the deformation of the surrounding soil is considered in the above-derived solution. On the other hand, the columnar inclusion also has deformation, which is related to the deformation of the surrounding soil. For the columnar inclusion, a unit cell is selected and the free body diagram for an infinitesimal element is shown in Figure 5. From this figure, the differential equation of columnar inclusion can be written as follows:  
\[\left\{