## 2 Numerical Analysis
### Determination of the Relieved Displacements
The stress distribution along the depth is initially sampled at the depths corresponding to the nodes for which the weight functions were evaluated. Figure 6: Example of the displacements induced by the tangential load component \(Q\), along the radial coordinate and for different positions of the Dirac delta function: (**a**) radial displacements and (**b**) tangential displacements. A piecewise linear stress distribution is assumed, and this can be found in agreement with the evaluation of the IFs at just the \(z\) points at which the Dirac delta function was iteratively applied. Thus, each of the IF integrals can be solved as shown below. \[ \begin{split} U_{*}&=\int_{0}^{H}\sigma(z)f_{U*}(z) \mathrm{d}z=\sum_{i}^{n}\int_{z_{i-1}}^{z_{i}}(\alpha_{1}+\beta_{1}z)(\alpha_{2 }+\beta_{2}z)\mathrm{d}z\\ &=\sum_{i=1}^{n}\alpha_{1}\alpha_{2}\Delta z+(\alpha_{1}\beta_{2 }+\alpha_{2}\beta_{1})\frac{\Delta z^{2}}{2}+\beta_{1}\beta_{2}\frac{\Delta z^ {3}}{3}\end{split}  \tag{4}\]  
This approach is already quite accurate due to the high resolution of the mesh, which implies a high density of the sampling points. However, an even more accurate calculation of the integral can be obtained by splitting the hundredth of a millimetre step if the stress distribution is introduced piecewise-linear with the interval points not coincident with the node positions. Therefore, the IFs were linearly interpolated on the same intervals, meaning that the calculation remains the summation of the trapezoidal product integrals (Equation (4)) except on an augmented set of intervals. The described approach for deducing the IFs directly from the FE simulations by means of the Dirac delta loadings allows to avoid any elaboration on the FE results themselves, with the exception of these simple linear interpolations, as shown in Figure 7b and discussed later. Only the final analytical integration of Equation (4) is required, which can be easily implemented in any coding language. According to the procedure explained here, the calculation could be possible only for specific values of the hole depth \(H\), namely those actually modelled with the FE simulations. In principle, this could may not be a significant limitation due to the high resolution of the FE model. However, given any \(H/R\) ratio, it is possible to obtain the evaluation of the IFs at the Dirac delta loading points by interpolating the numerical data of two depths which are multiples of the discretization steps and are higher and lower than the specific \(H/R\) ratio, as graphically explained in Figure 8a. The last point of the IF evaluations (along the horizontal red line) corresponds to the Dirac delta load just applied at the bottom of the hole, i.e., at depth \(H\). The IFs in this configuration can either be obtained by extrapolating the last two points, i.e., by moving horizontally (meaning at a fixed \(H\)), or alternatively through interpolation by moving diagonally, as shown in the figure, meaning that the load is imposed at the bottom while varying the depth \(H\). This latter approach is more accurate,  
Figure 7: Integration of the product of a generic IF, times the corresponding stress, obtained as the summation of the products of trapezoidal distributions: (**a**) integration resolution according to the FE model, (**b**) augmented resolution for a generic piecewise linear stress distribution. as interpolation is in principle preferred over extrapolation. This interpolation procedure is not performed when the \(H/R\) ratio is exactly coincident with a multiple of the mesh step. This is actually only effective for nominal cases or approximated measures obtained from rounded dimensions, which can, however, provide important reference cases for validation. The resolution based on 0.01 mm/1.00 mm is used for the initial depths, while larger steps are introduced later to avoid a heavy amount of data while preserving numerical accuracy.