## 2 Analysis Method in This Study  
In this study, SPH method, which is one of the most popular mesh-free methods, is adopted in order to analyze local fracture phenomena of concrete which are difficult to simulate by FEM. The physical quantity of each particle is evaluated by the weighted average of the adjacent particles using the kernel function, as shown in Figure 2, and the basic equation of the SPH method is written as Equation (1).  
\[<f(x_{i})>=\int_{\Omega}f(x_{i})W(x_{i}-x_{j},h)\,dx_{j}, \tag{1}\]  
where \(W\) is the kernel function, \(h\) is the smoothing length, terms inside angel bracket are SPH approximations, \(\Omega\) is the integral area within the support domain, and \(x_{i}\), \(x_{j}\) are the position of particle \(i\) and \(j\).  
Figure 2: Influence domain of SPH particles (kernel function).  
The smoothing length \(h\) is very important in SPH analysis. If it is too large, the accuracy of the analysis would be affected due to the smoothed-out details and local properties; on the other hand, if it is too small, the number of particles within smoothing domain may not be enough to calculate accurate response of particles [25]. Thus, in this study, an adaptive smooth length \(h_{a}\) is utilized. Since SPH method is to calculate weighted average values within the smoothing length, an adaptive smoothing length is better for the calculation of local deformation and stress concentration. The adaptive smoothing length \(h_{a}\) is changed during the analysis with the density of the particle, which is widely used in SPH analyses, as shown in Equation (2).  
\[h_{a}=h_{0}\left(\frac{\rho_{0}}{\rho_{current}}\right)^{\frac{1}{3}}, \tag{2}\]  
where \(h_{0}\) is the initial smoothing length, \(\rho_{0}\) is the initial density of the particle, and \(\rho_{current}\) is the density of particle calculated by SPH process, shown in Equation (3).  
\[\rho_{current}=\sum_{j=1}^{N}m_{j}W_{ij}. \tag{3}\]  
By the above calculation, the smoothing length \(h_{a}\) is decreased when the influence domain is compressed with the increase of density.  
When adaptive smoothing length \(h_{a}\) is adopted in SPH analysis, to obey [PERSON]'s third law and conservation of momentum, the following changes are applied to the calculation of kernel function, shown in Equation (4).  
\[<f(x_{i})>=\int_{\Omega}f(x_{i})W\left(x_{i}-x_{j},\frac{h_{i}+h_{j}}{2} \right)dx_{j}. \tag{4}\]  
By the above calculation, the interaction force between pairing particles by SPH analysis maintains the same.  
If two particles are positioned with smaller distance than \(h_{0}\) at the beginning of the analysis, we defined that the two particles are \"SPH-linked\", and the SPH-link between them will break when the distance between them is larger than \(h_{a}\), or the particles enter the crushing state, which is explained in the next section. Damage in solids is generally unrecoverable, and solid material, such as concrete, is not able to bond together due to the decrease of distance. Thus, if two particles are not SPH-linked at the beginning of the analysis, or SPH-link breaks by concrete fracture, no SPH link will be generated, even if the distance between particles become smaller than \(h_{a}\).  
To calculate the bond stress between concrete and steel anchor bar, the normal direction of the contact surface \(\mathbf{n}\) should be defined at the beginning of the numerical analysis. When the steel anchor is being pulled out, displacement \(\mathbf{dr}\) would occur within the pair of neighboring concrete particle and steel particle, as shown in the Figure 3a. Using the normal direction and the vector of the two particles, the slip along the tangential direction \(\mathbf{t}\) can be calculated as in Equation (5). Finally, by using the relationship between bond stress \(\tau\) and slip distance \(|\mathbf{t}|\), the bond stress can be obtained, as shown in Figure 3b and Equation (6).  
\[\mathbf{t}=-(\mathbf{dr}-\mathbf{dr}\cdot\mathbf{n}), \tag{5}\] \[\tau_{ij}=f_{b}(|\mathbf{t}|). \tag{6}\]