Residual Stress Determination with the Hole-Drilling Method on FDM 3D-Printed Precurved Specimen through Digital Image Correlation

Featured Application: This study shows the validation of an effective methodology for determining the residual stresses of components obtained with a 3D-printed polymer such as PLA based on the DIC and the hole-drilling method. Abstract: The hole-drilling method (HDM) is a common technique used for the determination of residual stresses, especially for metal alloy components, though also for polymers. This technique is usually implemented with strain gages, though other methods for determining the fields of displacements are quite mature, such as the use of digital image correlation (DIC). In the present paper, this combined methodology is applied to a 3D-printed PLA precurved specimen that is flattened in order to impose a bending distribution which can be considered known with a reasonable accuracy. The back-calculated stress distribution is in agreement with the expected (imposed) bending stress, however, a converging iterative procedure for obtaining the solution is introduced and discussed in the paper.


Introduction
The hole-drilling method (HDM) is a well known technique with a strong tradition, which allows to evaluate the residual stresses at the surface and within a certain initial depth, usually on the order of 1 mm [1][2][3].The HDM is usually implemented with a strain gage rosette that measures the relaxed (or relieved) strains.HDM with strain gages can be performed either as a through-hole (such as for thin sheet metal specimens), or more often as an incremental blind hole.Incremental hole-drilling requires that the strain gage readings are collected and saved for each step of the milling tool, and then converted into a certain distribution of stresses.These methodologies (through-hole and incremental blind hole) are covered by the ASTM E837 standard [4,5], which provides a reference for the implementation of the method, including best practices and calculation instructions to obtain the calibration coefficients.The implementation of this measurement technique can, however, require attention to certain details that are disregarded by the standard, such as the plasticity effect [6][7][8], possible small thickness of the plate in which the residual stress needs to be measured [9], the possibility of a chamfer at the bottom of the hole [10], and imperfect local planarity of the specimen, such as when the residual stress measurement is performed at the external surface of a cylindrical shaft [11].Another source of uncertainty is the hole's eccentricity; however, this can be corrected provided that the eccentricity itself is accurately known [12][13][14][15].Several analyses have been published in past years on different HDM errors and categorizing the different sources of error [16][17][18][19].A general description of the ill-posed HDM inverse problem was recently provided by Beghini et al. [20].In light of this ill-posedness property, regularization approaches and techniques such as in the recent ASTM E837 standard (2020 and 2013 editions) [4,5] are discussed, which are based on the Tikhonov regularization [21].For example, this approach is followed and implemented in several studies by Smit et al. [22][23][24].However, other approaches are possible, such as simple polynomial, spline filtering, or other techniques such as Gaussian Process Regression [25].However, the role of the regularization must be balanced in terms of the degrees of freedom introduced, which depends both on the quality of the measured data and the assumptions imposed in reaching the solution [20].
Regularization approaches can be tested on a known solution, which is very seldom available when operating with residual stresses.Of course, independent measures can be carried out and comparisons provided; however, even in this way, sources of incompatibility can arise between the solutions under comparison.For this reason, an interesting approach is to use a bending bench for validation, in which a well-known and simple stress distribution is imposed, and then retrieved, to test the quality of the inversion and its regularization [26][27][28].
The HDM can alternatively be implemented by recording and analysing the displacements near the drilled hole instead of reading the strains at specific locations, as with the strain gage rosette.This is basically done by combining the hole-drilling with one of two techniques: electronic speckle pattern interferometry (ESPI) [29][30][31] or the digital image correlation (DIC) [31][32][33][34][35][36].Recently, dedicated equipments have been proposed to implement the HDM in combination with the DIC technique [37,38], confirming the promising opportunity of this approach.
The DIC method can, in fact, be extensively and efficiently applied to different monitoring methods, especially of components produced by any additive manufacturing technique [39].In particular, this full-filed optical observation can provide the relaxation displacements of the entire surface of a component after the removal of clampings, which can in turn be interpreted to evaluate the residual stresses, as recently shown by Boruah et al. [40].
In the present paper, the HDM combined with the DIC technique (HDM-DIC) is implemented by initially proposing a finite element (FE) model based on the influence functions to obtain the analytical tool needed for the inversion.In order to impose a simple and known stress, a precurved specimen is then introduced and designed on the basis of the material properties.In principle, this resembles the use of a validation bending bench, as a certain flexural load can be imposed and then retrieved as the procedure output.The specimen was manufactured by means of 3D printing using a polylactic acid (PLA) material.The advantage of producing the specimen with this technique is twofold: first, the precurved specimen can be manufactured easily, reducing the required time and cost; second, the HDM-DIC can be tested on a 3D printed polymer, which can be a natural field of application of this technique, as despite the lower stresses, the (local) relieved displacements tend to be higher than in metals due to the lower stiffness.
The CAD file of the proposed specimen is available as a Supplementary Electronic Attachment to this paper for the dissemination and promotion of this methodology.

Determination of the Relieved Displacements
The hole-drilling problem requires the determination of calibration coefficients after having imposed the analytical form of the solution.The determination of these coefficients is usually performed by running different load cases and collecting the displacement results.In order to solve this problem with high resolution and accuracy within a reasonable computation time, a plane axisymmetric harmonic model is strongly recommended.As shown in Figure 1a, the radial displacement due to the relief of a circular volume and an equibiaxial stress distribution can be seen as the difference between the hole without the stress at the boundary and the hole (before drilling) with the stress at the boundary.The difference of these two loadings is equivalent to an inverse stress at the hole bound-ary and with no far-field stress.A similar approach can be applied to the shear loading (see Figure 1b), in which the tangential traction needs to be applied to the hole boundary along with the radial one and both stresses show an angular variability according to either cos(2ϑ) or sin(2ϑ), still with no far-field boundary conditions.The radial and the tangential displacements follow the same angular dependency.This ϑ variability for the two loadings is in principle only valid in a half-space geometry, which is often a reasonable assumption.On the contrary, specific analyses are required on more complex geometries, requiring dedicated 3D modelling, for example with the cylindrical components in [11].In this paper, the half-space was assumed to be representative of the analysed geometry under investigation.
(2 ) Initial condition with residual Q stress Zero stress at the hole boundary after drilling Following the concept described in Figure 1, it is evident that the load reduces to a traction at the cylindrical hole without far-field stresses.For this reason, the geometry reduces to a cylinder in a half-space with the far boundaries not loaded.Moreover, as the angular dependency is either constant or related to sinusoidal functions in terms of 2ϑ, the FE analyses can be performed with axisymmetric harmonic elements, e.g., PLANE25 in Ansys.According to this modelling, only a section with a plane through the axis of the hole is required, without the need to introduce any cumbersome and heavy 3D analysis.
The FE model and its features are presented in Figure 2. The half-space is modelled as a large cylindrical geometry with respect to the hole.Otherwise, the outer boundary can play a role, as discussed in [41].As there are small features and large global model dimensions in the same model, a strong mesh size variation is required from the hole region to the periphery.As shown in the figure, the hole radius was set equal to 1 mm.The region near the hole was imposed with a mesh of 0.01 mm and the mesh size was increased with a factor of ×2 at each transition, which was repeated five times from the hole region to the periphery.Despite only the radius R = 1 mm was analysed, a rescale of the problem can be applied, as discussed below, allowing different sizes to be easily considered with almost equal degrees of accuracy.
According to the ASTM E837 standard [4,5], the residual stresses can be well approximated by piecewise-constant distributions.Under this assumption, and with the principle of superposition being valid (as the system is linear), the effect at the surface in terms of displacements can be obtained as the sum of the single stress contributions σ i along the depth.Each of these needs to be multiplied by its weight U * i for a general component of displacement; this weight is the displacement generated by a unitary σ i and applied individually (see Figure 3a).In principle, this can be considered satisfactorily, especially when small steps are imposed.However, a more versatile approach can be obtained by replacing the summation with an integral of any continuous stress distribution.Thus, for each depth position there is a function f U * (z), here called the influence function, which resembles the unitary weight U * i for an infinitesimal depth increment dz (see Figure 3b).Examples of the concept of the influence functions (sometimes called weight functions, such as in fracture mechanics) can be found in [42,43] in the context of the HDM with strain gages.In principle, given the integral definition of such functions, the value at a certain depth z 0 can be found as the integral with a Dirac delta function which is singular at z 0 , as shown in Figure 3c, in which the units are recalled.The finite element (FE) technique offers a direct way to apply a Dirac delta function, as the load can be concentrated in single points (which are the nodes), while the values of these functions can be alternatively inferred from their cumulative results, as proposed by Beghini et al. [42].However, this requires the introduction of a model, such as a polynomial function, for each of them.In the present application, the use of a polynomial function would be complicated compared to the classical HDM with strain gage rosette due to the number of IFs being high.In fact, each point of the surface mesh grid starting from the hole rim up to a certain distance requires its own IFs in order to cover a large enough region around the hole.
The plane stress components at the surface of the specimen or the part to be investigated are σ x , σ y , τ xy .It is well known that for the HDM the use of alternative stress components is preferential, according to the following: Here, the P load is the equibiaxial stress, while Q and T are the shear stresses at 45 • and 0 • , respectively.Therefore, the Q and T loads are the same type of loading; for this reason, the FE simulation only requires P and one of the other two.Here, the Q load is reported in the initial Figure 1 and implemented in the following.The P load only generates a radial displacement, whereas the Q load implies both normal and shear tractions and in turn generates two components of displacements (the radial and the tangential ones; see Figure 4a and 4b, respectively).The use of point loadings according to the Dirac delta function can be easily performed for just the nodes along the depth; the displacements are recorded in terms of amplitude to be multiplied by either cos(0 × ϑ) = 1 for the P load or cos(2ϑ) or sin(2ϑ) for the radial and the tangential components related to the Q load.We used a resolution of 0.01 mm at the hole region for both the radial distribution of the IFs and the depth variable, as shown in Figure 4.In principle, the IFs are properly defined (i.e., not singular) for any point of the surface around the hole except for the rim point.When the Dirac delta function is applied at that point, the local displacement is singular.This problem is avoided by simply discarding the point at the very boundary of the hole.Thus, the first point, where the IFs are defined, is at r = 1.01.At this point, when z = 0 the Dirac delta function is close to but not exactly at that point; nevertheless, when the delta function is too close to the point under consideration, even when it is not superimposed, the solution is not expected to be extremely accurate, as the finite element introduces a local perturbation of the solution.This issue can be apparently circumvented by deducing the IFs from the cumulative displacement; however, this is inappropriate, as in an FE model any distribution of stresses can be seen as a superimposition of loads at the nodes, making it again a superimposition of the Dirac delta function solutions.An example of IFs along the radial path at the hole region is provided in Figure 5, in which the IFs are shown at single depth values of the Dirac delta function.From this graph, it is evident that large displacements can only be found close to the hole rim.For comparison with the strain gage method, the overall size of the usual grid is shown in the figure overlapped with the IF distributions.It is evident that the grid is relatively far from the hole rim where the IFs are lower.More precisely, the dependency of the displacements along the radial direction is not actually the IF itself, which on the contrary is the variability of these displacements depending on the depth position of the Dirac delta function.The other loading type (component Q) is shown in Figure 6, where for each point there is a radial displacement and a tangential displacement.The former follows cos(2ϑ), while the latter follows sin(2ϑ); thus, the maximum tangential displacement can be found at the angle ϑ = 45 • .The sign of the influence functions shown in Figures 6 and 7 are in agreement with a positive P residual stress component and a positive pure shear Q residual stress, while the scales of these distributions mainly depend on the Young's modulus and the depth of the hole.
Each IF actually depends on several terms: the material elastic properties (Young's modulus and Poisson's ratio, denoted as E and ν, respectively), the radial position at which the displacements are evaluated, the hole depth H, and finally the depth position z, which is the location at which the Dirac delta function can be applied to explicitly found the IF.The hole radius R = D/2 is another parameter that the IFs depend on; however, this latter geometrical variable is discussed later, and the hole depth is initially assumed to have R = 1 mm.Among these parameters, the only variability is the z coordinate, along which the integral needs to be performed after multiplying by the IF for the corresponding stress component according to the following.
When the stress directions are not the principal ones, a third stress component needs to be considered.This is the load T, which is again a pure shear, although now at an angular offset of 45 • with respect to the load Q.As this is again a shear stress, the same IFs need to be multiplied.These numerical integrations can be easily obtained according to the strategy reported in Figure 7a.The stress distribution along the depth is initially sampled at the depths corresponding to the nodes for which the weight functions were evaluated.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.
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, 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.The FE database was performed for a unique value of the Young's modulus E; this does not limit the generality of the problem, as all the displacements are linearly related to E. On the contrary, the variation of the Poisson's ratio modifies the numerical results of the displacements in a way that cannot be exactly predicted.For this reason, five Poisson's ratios were considered, and the full FE simulation set was repeated for these values.A piecewise linear interpolation was again performed for intermediate values of the Poisson's ratio (see Figure 8b).In order to avoid extrapolations, this procedure is not recommended for values of the Poisson's ratio outside the range of the simulated values, i.e., from ν = 0.25 to ν = 0.45, which is expected to include most materials of structural interest.
According to the procedure presented so far, the hole radius is R = 1 mm, and this would be a strong limitation if the radius were to be imposed.However, the scaling approach can be applied for this variable, which is an exact procedure.If the radius is different than unitary as modelled in the FE analysis, a simple scalar (dimensionless) parameter can be calculated for scaling the entire model: λ = D/(2 mm) = R/(1 mm), where D and R are the actual hole diameter and radius, respectively.
When the stress load distribution corresponding to a Dirac delta loading is applied, assuming that the problem is geometrically scaled by λ and considering the same applied traction, the problem simply scales with respect to the original configuration in which the diameter has a unitary radius.Accordingly, the corresponding points share the same solution in terms of strains; thus, the displacements are scaled by the same factor λ. The corresponding points are those obtained by scaling the coordinates, again by the λ factor.Reconsidering the IFs in this situation, the whole integral along the depth is higher or lower than the reference geometry by a factor of λ, as the z integration is scaled accordingly.This is due to the Dirac delta load being a limit distribution and the described scaling rule being valid for any instance of this limit.In conclusion, the displacements and the Dirac delta functions share the same scaling factor; therefore, the Dirac delta function with unitary integral generates the same displacements at the corresponding points, which in turn are just the influence function values.In the end, the IFs do not scale with λ; while the corresponding point coordinates scale instead.According to this rule, it is possible to retrieve the IFs at the surface points after scaling their radial coordinates from the original FE model.Finally, the integration of Equations ( 2) and ( 3) can be solved using the same method shown in Figure 7.
An example of the application of IF integration and determination of the displacements around the hole is shown in Figure 9, in which the parameters of the load, geometry, and material properties are similar to the experimental configuration reported later.

Inverse Determination of the Residual Stresses
The availability of the displacements at the upper surface when the residual stress is known can be referred to as the "direct problem".This tool actually allows us to solve the "inverse problem", which consists of finding the residual stress distribution that would produce the actually measured relaxed (or relieved) displacements.In principle, this turns out to be a minimization problem, which consists of finding the best distribution to obtain the minimum differences between the calculated and measured displacements.This minimization can be naturally solved according to the least-squares method due to the linearity between the residual stresses to be found and the relaxed displacements.As schematically shown in Figure 10, the displacements can be interpreted as a superposition of different stress contributions.In particular, if a piecewise-linear distribution of the stresses is assumed, then each step can be seen as the sum of two triangular stress distributions, which in turn can be seen as the scaling of a triangular distribution with the unity value at the maximum and the scaling factor equal to the stress at that vertex.
At each depth, the simple σ i shown in the figure is actually the term for stresses σ x , σ y , τ xy , or equivalently for stresses P, Q, T; thus, the following linear system can be imposed. [A] The matrix [A] in Equation ( 5) is 2N × 3n, in which n is the number of vertexes of the linear piecewise distribution, i.e., the number of stress steps +1, and N is the number of points detected by the DIC algorithm, for which two displacement components are available.As is obvious from the experimental side, N needs to be much larger than n such that [A] is a vertical matrix, otherwise the experimental data would be inappropriate with respect to the number of stress variables to be found.From Equation ( 5), each single column of the matrix [A] can be found by imposing a unitary component among P, Q, T and zero all the others, then repeating this for all three components of all the n stresses.For example, the first column is obtained by calculating the radial and the tangential displacements at the DIC points by imposing P 1 = 1 MPa, Q 1 = 0, T 1 = 0, P 2 = 0, . . ., T n = 0, then the second column is obtained by imposing Q 1 = 1 MPa and zero all the others, and so on.Each of these single distributions of stresses can be solved as a generic imposed stress by applying the algorithm presented above, in which the load is linearly varying and the step depths can be completely unrelated to the mesh points.Since the number of imposed equations (the measured displacements at the DIC points) is larger than the number of unknown stresses to be found, potentially even much larger, the system of equations is inconsistent (or even impossible).The natural way to find a solution is by imposing least-squares minimization, which can be obtained by calculating the pseudo-inverse of the matrix [A] and left-multiplying that by the column of the measured displacements.
Displacements (both rad.and tang.) Figure 10.Superposition of triangular unitary stresses, each scaled by a stress term which is the value at the vertex of the piecewise distribution to match the actual (measured) displacements at the upper surface.

Tensile Specimen for Elastic Properties
In principle, the proposed methodology is valid for determining the residual stresses on a specimen of any material provided that the constitutive relation remains within the linearity during the relaxation, i.e., when the hole is drilled.In this study, the material investigated is polylactic acid (PLA), which is a well-known and very common material in 3D printing, especially when combined with the fused deposition modelling (FDM) process.
A tensile test was initially performed to characterize the material, i.e., to find both the elastic properties (Young's modulus and Poisson's ratio) and the maximum tensile stress.As is evident in Figure 11a, the tensile specimen was 3D printed according to the plate type with a flat section, as proposed in the ASTM E8 tensile test standard [44].However, the cross-section had to be quite large in order to allow the application of both the extensometer and additional strain gages.The Young's modulus was determined by following the approach proposed by the dedicated ASTM 111 standard [45] and considering the extensometer measurement.Figure 11b shows the tensile curve up to large strains, where most of the section was fractured, as this fracture happened gradually.The very initial part of the curve was obviously of more interest for the present application.As is evident in Figure 11c, the first part of the curve was a ramp with a constant slope at least up to a certain limit, which was approximately at 20 MPa.After that stress, the curve shows a progressive curvature and the maximum stress of 35 MPa is reached.After selecting a section of the curve below the proportional limit, the Young's modulus of the printed material was obtained.This value was considered in the subsequent analysis instead of data from the literature, as this parameter is somewhat affected by the direction of the printing texture, the size of the filament, and other details of the printing process.As discussed before, the numerical procedure assumes a linear isotropic material, which requires the Young's modulus and the Poisson's ratio.Despite small variations of the latter, the material parameters produce little effect in terms of the back-calculated stresses; thus, an experimental determination of the Poisson's ratio was pursued.Two strain gage grids were aligned with the longitudinal direction and the transversal directions.These two grids were obtained with a strain gage rosette, in which only the two orthogonal strain gages were wired and measured (see Figure 12a).In principle, the longitudinal strain gage and the extensometer should provide identical measurements.As shown in Figure 12b, these two values were indeed quite similar, especially within the first strain range from 0 to 500 µε.After that, some small bending effects or any other source of misalignment could cause deviation between the two measures.The ratio between the transversal (absolute value) and the longitudinal strain slopes was calculated for the initial small strains, obtaining the Poisson's ratio.As shown in Figures 11 and 12

Precurved Specimen
The versatility of 3D printing was exploited in this study to design and manufacture a specimen allowing for the generation of a distribution of bending stress which could be considered known with sufficient accuracy to be used for a validation.The shape proposed here is a precurved plate which is flattened when three-point bending is applied (see Figure 13).This idea is not new, as a similar approach was proposed by Baldi [35] and previously by Schajer [46].However, their specimens were initially flat, and a curved seat was prepared to impose a certain degree of bending.
As our specimen is precurved, the seat of the specimen is flat and rigid enough to keep the specimen itself flat.The stiffness requirement can be easily fulfilled due to the printed material being quite flexible; thus, a simple planar support is enough to impose the bending of the specimen and ensure that it is ready for the HDM operation.
The design of the precurved specimen is simply obtained by imposing the opposite of the displacement to be experienced under the load.The beam theory under bending was assumed, as the width was intentionally quite small with respect to the length.Despite the material being quite flexible, the large displacement theory was not required, as the maximum imposed stress was limited up to the proportional limit (see Figure 11c).The midsurface sampling points shown in Figure 13a were obtained according to the following equation, which takes into account the anticlastic curvature according to the beam theory (y dependency): where z is the vertical upwards direction.A further vertical axis z ′ can be added along the depth direction with the opposite orientation (see Figure 13b).After this, the depth coordinate in the previous section should in principle be replaced with z ′ .The maximum height z obtained at x = ±L is the same as the maximum deflection of the cantilever beam formula: When bending is imposed and the entire z max is recovered at both sides, the bending load at the central line (x = 0) is where for the simple rectangular section I = ( 1 /12)bh 3 and W = ( 1 /6)bh 2 , with h being the specimen thickness (Figure 13b).More precisely, this latter bending formula (Equation ( 8)) was initially inverted to find the force F that generates a stress not exceeding the proportional limit, as discussed above.Equation ( 6) was used to find the points of the midsurface (Figure 13a).The surface generated with these points was then extruded in a CAD environment (Figure 13b) and the slots for the two screws were added.The specimen was printed with FDM (Figure 13c,d), for which a Tenlog Hands 2 Pro 3D printer equipped with a 0.4 mm nozzle was used, setting the extruder temperature to 205 • C and the layer height to 0.2 mm.The same printing parameters were used for the previously introduced tensile test specimen, and the longitudinal orientation of the precurved specimen, namely, the x axis, was the same as the tensile test direction to ensure the most similar material properties between the two specimens and that the loading was along the same printing direction.
As is evident in Figure 14a,b, and as introduced before, the precurved specimen was set to planar by two screws, generating three-point bending.The exact distribution of the contact pressure at the central point where the force 2F is applied depends on a number of factors, including the accuracy of the printing, the stiffness of the specimen seat, and the planarity of the seat.Controlling this distribution accurately is quite difficult; however, as the stress of interest is at the upper side, it is not as affected by the contact at the other side, as is evident in Figure 14c.This is investigated further below.

Stress Distribution
The stress distributions along the thickness of the plate specimen are in principle available from the beam theory.However, the distribution at the centre of the precurved plate is not perfectly in agreement with the theory, mainly because there is the contact with the load 2F in that region at the bottom surface, which implies a nonlinear distribution of the bending stress along the thickness (see Figure 15a).Despite this nonlinearity being more pronounced at the bottom side, it is advisable to apply the load and evaluate the stresses, not just the central point, but at a slight distance in order to avoid the effect of the contact and obtain a distribution more similar to that of the bending of a beam.Simple FE simulations allowed for estimating the position at which the solution is linear; for better accuracy, the stress distribution obtained with the FE model is preferred over the beam solution.Moreover, this side positioning, as shown Figure 15a, allowed for exploiting specimens with two distinguished measuring points.
More precisely, the contact at the bottom surface could be distributed in a different way than how it is modelled here.If the actual Poisson's ratio were (slightly) higher than the value determined above (Figure 12a), then the lateral sides would lose contact earlier than when the specimen is set flat, which would concentrate the force at the centre, as modelled here.On the contrary, if ν were lower than the measured value, the contact would remain distributed at the bottom middle line of the plate specimen.In this situation, the effect of Figure 15a is smeared along this length and locally reduced.It is worth noting that only the very initial distribution of the upper surface is of interest, as the hole depth is 1 mm or less, as shown below, and within this depth these perturbation effects are negligible.

Theoretical Framework and Implementation of Digital Image Correlation
Digital image correlation is a well-established technique to precisely measure displacements through image processing.The method is based on the presence of random features in the images, which can be tracked through a dedicated algorithm.These features can be naturally present on the specimen (e.g., surface roughness or texture) or can be added in the form of random speckle patterns to improve robustness.Due to the printing process, the 3D-printed specimen used in this paper presents a repetitive texture which is not suitable for DIC analysis; thus, black speckles were sprayed on the top white surface before the measurements.When such a random pattern is present, it is possible to perform DIC measurements by comparing the deformed image to the undeformed one.Many different DIC algorithms are available in the literature, most of them based on a two-step procedure involving rough displacement estimation at the pixel level through correlation followed by fine subpixel measurements through interpolation schemes [47].In this paper, the Digital Image Correlation and Tracking algorithm was used [48].To accomplish the first step, a measurement grid can be arbitrarily defined on the image.An example grid is shown in Figure 16a with red dots.Then, a square subset is defined around each measurement point, e.g., the red box in Figure 16a.The algorithm then iteratively moves the subset in the deformed image, as exemplified in Figure 16b with the green box.The correlation between the moved subset and the original subset in the reference image is then computed using the following equation [48]: where u and v refer to the iterative tentative displacements along the horizontal and vertical directions, respectively.The indexes i and j refer to the image pixels, while r and d represent reference and deformed image subsets, respectively.Additionally, the terms r and d represent the mean value of both reference and deformed subsets, and are defined for normalization purposes.The correlation result C is a matrix that contains information about the correlation value, i.e., a scalar between 0 and 1, for each tentative displacement u,v.It is worth noting that because the images are discrete matrices, the values of u and v at this step must be integers to allow for the selection of the pixel intensity values.The resulting correlation map can be used to detect the maximum correlation at the pixel level, which represents the rough displacement between the reference and deformed images.To improve the method's sensitivity, it is possible to refine the measurements for subpixel accuracy.In particular, a sub-array of C is defined by selecting eight points around the rough maximum, for a total of nine measurement points.These points are then used to define a two-variable quadratic polynomial fitting function, which allows for analytical detection of its maximum.The peak of this polynomial function is assumed to be the actual maximum of the correlation map, and the fitting procedure allows for an accuracy between 1/10 and 1/100 of a pixel.

Experimental Setup
The experimental setup was implemented on a milling machine to ensure easy and well controlled translational displacements (see Figure 17).The specimen was flattened on a plate, which was fixed to the (movable) base of the machine.The DIC camera was fixed to the upper frame of the machine so that it was fixed along any direction, and a light ring was positioned at the extremity of the camera lens to improve the local lightening and reduce the exposure time.A global shutter DMK 38UX540 camera (The Imaging Source) was used, having a resolution of 5320 × 4600 (24.5 MP).The field of view was set to be approximately 8 × 7 mm 2 , for a scaling factor of about 1.549 µm/pxl.The vertical distance between the camera lens and the specimen was initially fixed to ensure the best focus at the hole's upper surface.The milling tool was initially brought into contact with the upper surface of the specimen.In order to avoid slight damage to the surface before the actual cutting, a calibrated thickness was placed between the milling tool and the specimen surface, while rotating the former for detection.This thickness was taken into account to finally set the column vertical position corresponding to the zero depth of the milling tool.This latter was then set to contact while not rotating and a depth indicator was introduced at this position, which was fixed to the vertically sliding column of the machine, leaving it fixed with respect to the milling tool's vertical position.At this status the depth indicator was zeroed, as represented in Figure 17.Despite the tip of the depth indicator being far from the milling tool, the rigidity of the machine frames ensured the same translational displacement between the indicator and the milling tool.After this zeroing, the cutting tool was initially brought up (along with the depth indicator), then the base was moved leftwards by a length equal to the distance between the camera and the milling tool axes.After taking the initial DIC picture without the hole, the base was again repositioned below the milling tool.The drilling operation was then started with accurate control of the drilled depth due to the depth indicator.The rotating speed of the milling tool was imposed at 200 rpm and the feed rate was 0.1 mm/s.These values were kept low to avoid local overheating and melting of the specimen material.After drilling, the base was again moved to the position that the specimen was below the camera, and the final DIC image (with the hole) could be obtained.Despite the depth being measured with a good accuracy, the diameter of the generated hole is only known as a nominal value.An accurate acquisition of the diameter is obviously required, as it is both the main parameter of the HDM and the scaling length of the DIC displacement measurement.After the drilling and the DIC acquisition (Figure 18a), an optical device was used to observe the drilled hole with a thin frame (Figure 18b).This latter can be translated and the frame movement was associated with a dial gage with a resolution of one hundredth of a millimetre (Figure 18c).This device was available with the RESTAN instrument supplied by SINT Technology Srl, Florence, Italy.By moving the frame with respect to the specimen, it is possible to align one of the two lines (horizontal or the vertical) with the lateral rim of the hole, and then the gage translation to reach the other side can be easily measured.This can be repeated for the other line, allowing the average between the two diametral lengths to be finally calculated.The produced holes showed accurate roundness, confirming the correct applicability of the cutting parameters to the PLA.

Results
This HDM-DIC procedure was applied to different holes manufactured with an offset of a few millimetres with respect to the axis of the specimen, as shown in Figure 18a and motivated by Figure 15.In this analysis, the results of two tests at two different nominal depths of 1 mm and 0.5 mm are presented.The resulting hole diameters were 1.805 mm and 1.815 mm, respectively, and were obtained with the typical milling tools used for the HDM combined with the strain gages.The results obtained here can be directly compared to the stress distributions shown in Figure 15b at the corresponding initial depths, either 1 mm or 0.5 mm for the two cases.
In principle, the stresses induced in the present specimen should not be defined as "residual", as the specimen is loaded remotely instead of actually carrying those stresses.More precisely, the initial depth could be affected by both residual stresses induced by the manufacturing operations and by the imposed bending stress.The possible residual stresses were not evaluated in this study; they were assumed to be relatively low with respect to the bending, as the spring back was almost imperceptible after manufacturing.
Taking into account the shape of the IFs presented before, it is quite evident and obvious that the largest displacements are close to the hole rim and rapidly vanish.For this reason, not all of the observed points with displacements measured by the DIC should have the same significance when the stress inversion is performed; in particular, the points quite far from the hole should be excluded from the evaluation.In order to apply this, a simple approach could be to define an annular region of selection.However, according to the solution, the displacements could be high along some directions, while they could be rapidly falling along others.The strategy proposed here is an iterative procedure in which the remote points or those with little displacements are iteratively disregarded.Initially, all of the points are considered and a first solution is found.According to this first estimate of the stress distribution, the displacements are calculated at the DIC points themselves, then a selection criterion is implemented.A (dimensionless) accuracy factor f is initially found for each DIC point: where U r,N , U ϑ,N , U r,E , U ϑ,E are the radial and tangential displacement components, respectively obtained numerically (N), from the first solution, and the experimental ones (E).
According to this, a first selection of the points was taken on the basis of this accuracy factor.An effective choice was 0.667 ≤ f ≤ 1.5.This criterion was combined with the selection of those DIC points having U N ≥ 0.05 U N,max , where U N,max is the maximum displacement norm (numerically evaluated) of the whole observed field.The application of this criterion provides a deselection of both the outliers, which can have a significant weight (especially when the measured displacements are very different than expected), and of those points which are too far from the hole rim, thus showing small displacements.It is worth noting that the far points may have large weights, as the area of the peripheral region is significant with respect to the area close to the hole region.Moreover, some spurious effects (such as nearby constraints or lateral edges) can introduce inaccuracies between the measure and the model, even if at very small scale of the displacements, especially far from the hole.
After the first selection of the points, the solution was calculated again and the selection was repeated.This new selection was performed starting back from all the initial points, as some points which were initially disregarded could turn out to be well correlated as the solution evolves.This iteration was repeated and the procedure was found to converge quickly to a stable selection of points in approximately 10 iterations.
The results of the first hole (nominal depth 1 mm) are presented in Figure 19.The selection of the points at the end of the iterations is reported in Figure 19a, while the iterative solution is shown in Figure 19b; the same graph also presents the expected stresses for validation.The displacements at the end of the iteration procedure are shown in Figure 20, which presents both radial and tangential components, along with both the numerical and the experimental distributions, which are the variables involved in the evaluation of the accuracy factor in Equation (10).The same analysis was repeated for the second hole with nominal depth 0.5 mm.In this situation, the solution was less stable in terms of component gradients because of the smaller depth.
For this reason, the same procedure was equally repeated after imposing the condition of a uniform stress for each component along the entire depth of the solution.Recalling Figure 10, this was easily implemented by imposing the same levels of the superimposed triangular distributions; in other words, the same stress was shared by the two unitary distributions covering the entire depth, for each of the P, Q, T components, leading to just a single unknown for each of them.
The obtained stress distribution and selection of the points are shown in Figure 21.In this situation, the procedure converged even faster (approximately 3 to 5 iterations were enough to find a stable result), and again the obtained stresses were found to be quite in agreement with the expected ones.
The components of relieved displacements are shown in Figure 22.The pattern of the selected points changes slightly, and because of the smaller depth, the magnitude of the displacements is somewhat lower than the example.However, these displacements were still well detectable with the present settings by means of the DIC technique even for this shallower hole with 0.5 mm depth.The same approach could be proposed for the previous test with higher hole depth, as the slopes of the obtained stress distributions are not in perfect agreement with expectations.A similar stability and solution convergence rate was obtained for the 1 mm hole depth (not shown for brevity).However, the uniform stress combined with the proposed DIC point selection procedure would be more interesting for shallow hole depths even lower than 0.5 mm.In fact, an accurate evaluation of the stresses could be more useful quite close to the upper surface, where the residual stresses may be higher and most damaging, especially considering that the imposed uniform stress is more acceptable on a smaller depth.
Other techniques can be applied to smooth the DIC relived displacements, as proposed in Refs.[31,35].On the contrary, the method proposed in this study uses just the raw displacements to avoid any artifacts due to displacement smoothing or further data postprocessing.However, a certain number of acquired points need to be disregarded from the initial whole acquisition, according to the iterative procedure, to have a stable solution which was found to be accurately matching with the expected stresses.

Conclusions
The present paper proposes a high-resolution numerical method for calculating the relieved strains based on axisymmetric harmonic finite elements and direct determination of the influence functions (IFs).This allows for an efficient solution of the hole-drilling method (HDM) combined with the DIC technique.In addition, this paper proposes a precurved specimen which can be easily 3D printed with any additive manufacturing technique, e.g.fused deposition modelling (FDM).This specimen can be flattened on a rigid ground and the induced stresses can be easily determined with the beam theory, or more accurately with a simple FE analysis to consider the small variations with respect to the bending formula.In the experimental part of this study, specimens were printed with polylactic acid (PLA), with the FDM technique, and evaluation of the induced bending stresses was performed on two holes of different depths.The results were found to be in good agreement with the expected bending stresses; however, a deselection of the measured DIC points was required, and thus an iterative procedure was proposed which turned out to be efficient and quickly converging.An interesting result was obtained by imposing a uniform stress distribution for a hole with lower depth.The procedure converged even faster and the final distribution was in good agreement with the imposed stresses, confirming that this technique can be applied to this kind of material even at shallow depths.

Figure 1 .
Figure 1.Relieved displacement distributions obtained as the difference between the holed and the non-holed configurations for (a) equibiaxial stress and (b) shear stress.

Figure 2 .
Figure 2. Details of the axisymmetric harmonic FE model: (a) identification of the modelled plane and definition of the hole radius, (b) overall view of the FE model and dimensions, (c) number of nodes and elements, and element sizes, (d) mesh transition method.

Figure 3 .
Figure 3. (a) Piecewise-constant distribution of the residual stress and displacement obtained as a summation of unitary weights.(b) Integral definition of a generic influence function.(c) Extraction of the influence function values by introducing a Dirac delta function.

Figure 4 .
Figure 4. Dirac delta loadings at different points along the depth for (a) P and (b) Q components and numerically defined influence functions at the surface nodal points.

E
and slightly on ν

Figure 5 .
Figure 5. Example of the displacements induced by the load component P along the radial coordinate for different depths of the Dirac delta function, and position and size comparison with a rosette strain gage grid.

Figure 6 .
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.

Figure 7 .
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.
Piecewise linear interpolation of the databases on the Poisson's ratio

Figure 8 .
Figure 8.(a) Definition of the simulated configurations in terms of hole depth and Dirac delta function applications points along with possible intermediate interpolations; (b) different databases with the same resolution and loading points obtained with different values of the Poisson's ratio.

Figure 9 .
Figure 9. Example of applying the procedure to obtain the displacements at the top surface: (a) stress distribution assumed as known, (b) definitions of the stress components and orientation of the output displacements, (c) radial and (d) tangential obtained displacements.

Figure 11 .
Figure 11.(a) Tensile test specimen in PLA with both extensometer and strain gages.(b) Tensile curve up to almost complete fracture, obtained with the extensometer measurement.(c) First part of the tensile curve, showing the Young's modulus determination and the proportional limit.

Figure 12 .
Figure 12.Acquisition of transverse and longitudinal grids for the determination of the Poisson's ratio: (a) wiring of two orthogonal grids of a strain gage rosette, (b) ratio between the different slopes with respect to the extensometer.

Figure 13 .
Figure 13.Design and manufacturing of the precurved specimen: (a) definition of the plate midsurface by sampling points, (b) 3D solid obtained in CAD environment, (c) 3D printing operation with small supports to allow the curved shape, and (d) completed printing operation, showing the supports at the bottom surface to be removed.

Figure 14 .
Figure 14.Bending loading on the precurved specimens: (a) initial positioning of the specimen and the screws, (b) application of the vertical forces generated by of the lateral screws and force reaction at the centre, and (c) bending distribution simulated by FE and verification of the planarity of the shape under loading.

Figure 15 .
Figure 15.Stress distributions along the thickness of the plate: (a) at the centre of the plate, showing the perturbation effect due to the bottom contact and (b) at a slightly lateral position, where the solution is more similar to the beam theory distribution and is expected to be more reliable.

Figure 16 .
Figure 16.DIC procedure, definitions of reference (a) and shifted (b) subsets, and introduction of the displacement components with respect to the hole centre coordinate system.

Figure 17 .
Figure 17.Picture of the experimental setup, including the milling tool, DIC camera, and depth indicator.The movable base allows the specimen to be positioned either below the DIC camera or below the cutting tool axis.

Figure 18 .
Figure 18.Measuring the hole diameter: (a) holes at the upper surface of the specimen, (b) observation of the thin movable frame, and (c) dial gage with a resolution of one hundredth of a millimetre.

Figure 19 .
Figure 19.Residual stress determination for the hole with depth 1 mm: (a) selection of the points with accurate correlation and large displacement, (b) obtained stresses at the different iterations (lighter lines) and at convergence (solid lines), and validation with the expected stresses (dashed lines).

Figure 21 .
Figure 21.Residual stress determination for the hole with depth 0.5 mm: (a) selection of the points with accurate correlation and large displacement and (b) obtained stresses at the different iterations (lighter lines) and at convergence (solid lines), and validation with the expected stresses (dashed lines).