Stress Concentrations in Skew Pressurized Holes: A Numerical Analysis

This work provides a numerical analysis of the stress concentration factor in an elastic solid containing non-aligned and non-concurrent circular holes subjected to an internal pressure. Using Finite Element Analysis, a variety of systems with a range of geometric configurations and loading conditions were simulated and the trends observed were analysed using a qualitative and statistical analysis in order to determine the correlation between these factors and the maximum stress concentration in these systems. Furthermore, simple empirical models were calibrated on the simulation results and used to plot reference graphs which may be employed to predict the stress concentration factor of these systems according to several geometric parameters and loading conditions. The numerical results and the empirical models presented here also show good agreement with previously derived analytical results based on 2D models and are expected to provide a useful tool for the designer of such systems especially for fatigue problems in fluid-power systems, where a straightforward evaluation of the stress concentration factor in pressurized hole systems is needed.


Introduction
In this work, an analysis of an elastic solid possessing pressurized skew (non-concurrent and misaligned) circular holes was conducted.The maximum elastic stress is retrieved for a number of geometric configurations, including varying proximity, size and orientation of holes with respect to one another (i.e.angular rotation), and load-cases used to estimate stress concentration factors useful for the design of many industrial components.In the past, the case of a structural component containing a pressurized hole close to another or to a boundary and the resultant stresses derived from the interactions between them has been the subject of several studies due to the presence of these configurations in fluid power systems.The pressure imparted by the fluid could lead to material failure and fatigue cracks which originate and propagate from regions of high stress concentration such as those found around the periphery of holes and thus an analysis of the stresses found in these systems is essential in order to design and assess their structural integrity.
Several relevant studies are collected in Peterson's Stress Concentration Factors [1], although they are mainly concerning concurrent perpendicularly-aligned pressurized holes in thick tubes of various diameters, by far the most frequently encountered type of system and these could be adapted to perpendicular concurrent holes of various diameters.Moreover, other analytical approaches have been used to examine this situation in various forms and configurations; the majority of which involve a 2D elastic analysis.One of earliest considered cases involved a pressurized eccentric hole in a circular disc [2][3][4][5][6][7][8][9].Using bipolar coordinates, analytical solutions have been proposed for the calculation of the principal stresses found in such systems, which are expressed in terms of the geometric parameters and level of applied pressure.
Another commonly studied configuration is a circular hole near a straight-line boundary of a half-plane [2,3,10].Both these forms may be used to analyse the stress concentrations of a hole close to the periphery of a block of material.
Studies have also been conducted to investigate the interactions between two circular holes in close proximity to each other in an infinite or quasi-infinite plane [11,12,[21][22][23][24][13][14][15][16][17][18][19][20].The effects have been studied for cases involving equally and variably pressurized holes as well as holes subjected to uniaxial and biaxial stresses.Other related cases which have been considered in literature [25], involve the effect of hole shape (in cases of elliptical [26] or polygonal holes [27][28][29]) on the stress distribution mapping of the system, the influence of reinforcement on the periphery of the hole [30], periodically-replicated hole arrangements [31,32] and coalesced holes [33], as well as holes subjected to thermoelastic [34,35] and external stress loading [36,37].These problems have been tackled using a variety of analytical and numerical methods which have yielded both exact and approximate solutions, with varying degrees of complexity and mathematical robustness.Scirè Mammano and Dragoni [16] attempted to condense and simplify this problem by approaching it from the perspective of a bipolar coordinate representation, expressed in terms of two normalized geometric parameters; the relative diameters of the two holes (which tends towards infinity in the case of a straight-line boundary) and the ligament thickness (i.e.shortest distance between holes) to hole diameter ratio (see Figure 1a,b,c).This, in turn, allowed for the plotting of maximum stress vs geometric parameters graphs which provide a simple and direct reference to the engineer for the manufacture and design of such a system without necessarily delving into the mathematics and complex derivations used to obtain these plots.
The plots presented in this previous work [16] may only be used to predict the maximum stresses in 2D systems, i.e. systems where the holes are orthogonal to the same plane.However, while this is sufficient as a point of reference for the design of hole systems on a plane, it does not give a complete picture of the stress distributions and concentrations found in a system where the holes are un-aligned in 3D space.Such systems are commonly found in industrial settings, with the most typical example being that of the hydraulic manifold (see Figure 1d).
To estimate the fatigue strength of these systems, it is imperative to consider the problem from a 3D perspective, rather than from a 2D viewpoint; a situation which has not been covered to the same extent as the 2D configuration in literature.Indeed, relatively few works have examined the interactions of cylindrical hole systems in blocks of materials.Badr et al. [38,39] examined systems consisting of two orthogonally aligned cross-bores while other studies have considered the effect of single holes intersecting cylindrical and rectangular beams [40,41].
Other studies found in technical literature have investigated primarily the fatigue performance of three-dimensional hydraulic components [42][43][44][45][46][47].In view of this, in this work, we will tackle the problem presented by the placement of skew holes in close proximity with each other and the maximum stresses arising from the introduction of equal and unequal pressures within these cavities.More specifically, we investigate the system shown in Figure 2, which is the 3D equivalent of the planar configuration covered by Scire Mammano and Dragoni [16] (shown in Figure 1c), with the addition of a relative angle of rotation between the two cylindrical holes.Using Finite Element Analysis (FEA), the effect of misalignment of the holes on the maximum principal stresses of the system, coupled with the normalized geometric parameters of the holes, was analysed and used to obtain reference plots for these systems.

Geometry of the System
The simulations were conducted using the ANSYS16 Multiphysics platform.In order to obtain the internal stresses of a system corresponding to an infinite or very large block of material, i.e. study the system without boundary interactions, the system was modelled as a large cylinder with two cylindrical holes bored right through the central region, as shown in Figure 2a.The circular faces of the cylinder were aligned orthogonally to the y-axis, while the first circular hole, with diameter, d0, was cut along the z-axis.The second hole, with diameter d1, was cut at a variable angle, θ, with respect to the first hole, as shown in Figure 2b.The shortest distance between the periphery of the holes, or the 'ligament' thickness may be defined by the variable t, which results in a distance between the hole centres of d0/2 + t + d1/2 (see Figure 2c).The parameter d0 was kept constant at 20 mm for each simulation while the parameters d1 was varied: 5, 10, 15 and 20 mm, and the parameter t: 2.5, 5, 10 and 20 mm, resulting in normalized d1/d0 ratios of 0.25, 0.5, 0.75 and 1, and for t/d0: 0.125, 0.25, 0.5 and 1.This range of parameters encompasses the main cases for adjacent hole systems typically found in industrial applications [16].The parameter θ was varied as follows: 0°, 15°, 30°, 45°, 60°, 75°, 90°, with the system having a θ value of 0° being equivalent to a cylinder with two aligned holes orthogonal to the same plane (same case as in [16]) and the structure having a θ value of 90° possessing two holes which are aligned perpendicularly with respect to each other.The dimensions of the overall cylindrical block of material were defined according to the hole and ligament dimensions.A diameter, D, of 260 mm, equivalent to 13(d0), was chosen, while the height, H, was set to 4(d0 + t + d1).These dimensions were chosen to achieve accurate results, thanks to convergence test not reported here for the sake of brevity, whilst maintaining a reasonable level of computational efficiency.

Simulation Method
The system was meshed using SOLID187 elements, a higher order 3D tetrahedral, 10-node element with quadratic displacement behaviour.In order to maximize computational efficiency and optimize the number of elements per simulation, the system was initial meshed using the smart sizing option available in the ANSYS software at a fine level and then followed up the additional refinement of the mesh at the periphery of the holes and the central segment of the ligament (intersection).This method typically resulted in a system containing a minimum of 1.2 million elements.The results of the mesh convergence analysis are presented in the Supplementary Information.As shown in Figure 3, the system was constrained from three points in space; at one of the bottom face edge points a displacement fix in the x-, y-and zdirections was applied, at the corresponding opposite point in the bottom cylinder face, a fix in the x-and y-directions and in the point directly above it on the upper cylinder face, a fix in the x-direction only.This constraint method was used in order to minimize the constraints within the system allowing it to deform freely whilst ensuring that rigid body motion is inhibited.The complete set of geometric parameters, material properties and loadcases used to model these structures are summarized in Table 1:

Results and Analysis
The maximum 1 st principal stress, σ1p,max, was extracted from the central region of each simulated structure, so as to eliminate possible boundary effects which are not of interest in this study.In order to validate the simulations results, the maximum 1 st principal stress values obtained from the FE analyses with θ = 0° were compared with those obtained by Scirè Mammano and Dragoni [16] both for analytical and FE simulations of 2D systems of holes in close proximity under plane-strain conditions (shown in Figure 1c).A comparison of these results is shown in Appendix 1, which shows very good agreement compared to the FE simulations in particular, indicating that the chosen size and shape of the material block and stress value extraction method, are appropriate to obtain values which are comparable to those one would expect from a system with two holes in an infinitely-thick block of material.For the analytical results, we found the same discrepancies observed by the previous authors [16], especially for systems with low t/d0 ratios, where it was reported that the analytical expression may overestimate the stress concentration factor by up to 20%.
The results obtained from the three sets of simulations are presented in this section, where the maximum 1 st principal stress is first normalized against the applied internal pressure p to obtain the stress concentration factor, Kt, ( ), and is then plotted against θ for the various set of dimensionless parameters t/d0 and d1/d0.For cases where pressure is applied on one hole only, the parameter p indicates the magnitude of the pressure applied on the pressurized hole only.In the following subsections the results obtained for the three loadcases are analysed and the effect of the variables considered in the FE analyses are discussed.In addition, the findings of an Analysis of Variance (ANOVA) conducted on the separate datasets are provided in the Supplementary Information, as well as empirical models obtained by curve fitting techniques to predict the stress concentration factor in these systems.The ANOVA highlights what are the variables that influence the most the stress concentration factors from a statistical standpoint.
This powerful technique, which is part of the Design of Experiments method, born mainly for experimental tests, [49] could be profitably used also in numerical analysis [50], in order to find out if all the variable considered in the problem have a significant effect on the response.

Loadcase A: Equal pressure applied on both holes (p0 = p1 = 1 MPa)
The results obtained from the FE simulations conducted using Loadcase A conditions are presented in Figure 4, while a diagram showing the normalized 1 st principal stress distribution found for the system with the following parameters d1/d0 = 0.5, t/d0 = 5 and θ = 30° subjected to the same loading conditions is presented in Figure 5.It is evident that the parameter θ has a strong effect on the level of maximum stress, with the level of stress being at its highest level in the cases where θ = 0° and lowest where θ = 90°.This effect is particularly pronounced for smaller values of t/d0 ratios, where in some cases the maximum stress observed in the θ = 0° is double that for the corresponding system with θ = 90°.Moreover, in accordance with findings of previous studies, the maximum stresses were attained when the t/d0 ratio is very low.It is also apparent that as the t/d0 ratio increases, the Kt value tends towards 1, which is the value obtained when a uniform pressure is applied on a single hole in infinite elastic space.This indicates that the interactions between the pair of holes become negligible as this ratio increases.Finally, it can also be observed that the d1/d0 ratio has a negligible effect on the maximum stress value obtained for this particular loadcase.This observation was also confirmed by the ANOVA results which showed that this variable has a negligible effect on Kt.The results of the numerical analysis can be used to produce reference plots which may be employed to predict the normalized maximum first principal stresses of similar 3D non-aligned pressurized holes under the same loading conditions.Using optimization curve-fitting techniques, a 3D empirical equation was fitted over the numerical results and used to correlate the variation of stress with the changes in geometric parameters t/d0 and θ for Loadcase A (d1/d0 was not considered since it was shown to have a negligible influence).Based on the results of the ANOVA and the trends observed from an individual analysis of the 2D plots of the geometric parameters against maximum stress, the following function was fitted over the results obtained for the cases where p0 = p1 = 1 MPa: Eq. 2 This equation may also be expressed approximately in the following simplified and compacted form with a quasi-equivalent R 2 value:

Eq. 3
The 2D curves which may be used to generate reference graphs for this particular loadcase based on this simplified empirical model, are presented in Figure 6.

Loadcase B: Pressure applied on larger hole (d0) only (p0 = 1 MPa, p1 = 0 MPa)
As one may observe from Figure 7, the results for applying pressure to one hole only are distinctly different from those obtained for pressurizing both holes.While the trends obtained when changing the parameter θ are similar, with the maximum stress level decreasing as θ increases, in this case, the d1/d0 ratio is no longer an insignificant geometric variable and has a pronounced effect on the determination of Kt.In fact, it is evident that as the d1/d0 ratio increases, the maximum stress level also increases.On the other hand, for the t/d0 ratio, the opposite trend is observed, with Kt increasing as this ratio decreases.Interestingly enough, the results shown in Figure 7 also indicate that in cases where the distance between the holes is equal to or greater than the radius of the largest hole, i.e. t ≥ d0/2, the effect of θ on the maximum normalized stress is extremely small.Moreover, all systems with t/d0 values greater or equal to i.e.where the maximum 1 st principal stress is equal to the applied uniform pressure, which can be taken as an indication that the distance between the holes in these cases is sufficiently large to ensure that there is almost no interaction between them.Similarly, to Loadcase A, an empirical model was fitted over the numerical data points and used to produce reference plots to predict the 1 st principal stress normalized over the applied pressure.As stated previously, the geometric parameter d1/d0 has a significant effect on the maximum stress level in this case and, therefore, it was included as an additional variable for the curve-fitting procedure.As a result of this, for Loadcase B, a 4D model of the form shown in Eq. 4 was fitted over the data points obtained from the numerical simulations: An R 2 value of 0.995 was obtained for the model shown in Eq. 5 and reference plots showing the model plotted along with the numerical data points are presented in Figure 8.

Loadcase C: Pressure applied on smaller hole (d1) only (p0 = 0 MPa, p1 = 1 MPa)
The trends observed for this loadcase are similar to those obtained for Loadcase B (see Figure 9).It is clearly evident that lower maximum stress values are obtained when only the smaller hole is subjected to uniform pressure, which was to be expected.Due to his factor, the influence of θ on Kt is also diminished in comparison to the other cases.Finally, same as in Loadcase B, all geometric configurations where t/d0 is greater or equal to 0.5, exhibited Kt values lower than 1.3, indicating that there is minimal interaction between the two holes.The following empirical expression was obtained (Eq.7), with an R 2 value of 0.996, and reference graphs are shown in Figure 10.The term b obtained from the curve fitting was extremely small and therefore can be eliminated from the final equation.

Discussion
It is envisaged that the empirical models presented in this work will be useful as reference guides for the design of skew pressurized hole systems such as those found in hydraulic manifold systems.While these models only provide an estimation of the actual maximum stress values present in these systems, they highlight very clearly the trends resulting from changes in the geometric parameters and demonstrate how one may vary these parameters in order to lower the stress concentrations within the system.Furthermore, it can also be observed from the plots in Figures 6, 8 and 10 that at low t/d0 values, i.e. the geometric conditions which result in the highest normalized maximum stress values, the empirical models mostly err on the side of overestimating the stress concentration factor with respect to the numerical results.
This indicates that these empirical models may be reliably used as a pre-design tool to estimate the safety factor for multi-hole systems.In this work, the stress concentration factor is obtained by dividing the 1 st principal maximum stress over the applied pressure, however, using the data obtained for the 2 nd and 3 rd principal stress values obtained for the same point (presented in the Supplementary Information), one can easily apply the findings of this work to determine the It is also important to highlight the fact that in all systems studied here, the maximum stress value was to be found in the 'ligament' region between holes.This was to be expected, of course, however it is worth mentioning that while in the majority of cases, the maximum value was to be found at a single point in the ligament region along the y-axis at the exact centre of the cylindrical block of material (or extremely close to it), in some cases the stress concentration factor was found at a point significantly displaced from the centre or in some cases, where θ was equal to zero, at two points.The latter situation has also been observed and predicted by analytical models for 2D systems in previous studies found in literature [16][17].
Considering the results reported in Figures 6, which is the most common loadcase, it is worth noting that the diameter ratio does not play an important role, which means that the designer should focus mainly on optimizing the shortest distance between the holes, which should be increased in order to lower the stress concentration.The exponential-like curve shows as well that a small increase in t/d0 results in significant stress concentration reduction and, as expected, the best-case scenario occurs when the holes/ducts are perpendicular to each other.Focusing on Loadcase B and analysing the results in Figure 8 it is relevant to report that, for every parameter condition, when t/d0 exceeds 0.5 the other variables become irrelevant and the maximum stress concentration becomes less than 1.2.This conclusion is valid to the same extent for Loadcase C (see Figure 10) which means that irrespective of other geometrical parameters, when only one hole is pressurized, it is sufficient to have the shortest distance between the holes larger than the radius of the larger hole and the maximum stress value becomes equal to the applied pressure, which could become a practical rule of thumb that helps the designer.Loadcase C highlights as well also a lesser sensitivity to the angle for systems where the diameter ratio is lower than 0.5, since in this case the curves almost collapse into a single line and therefore the angle become unimportant in determining the maximum stress level.Moreover, for Loadcase B and C when the ligament becomes comparable with the larger diameter (t/d0~1) the stress concentration factor becomes unitary meaning that the two holes are not interacting anymore.
Finally, it is also worth mentioning that the reference plots and models presented in this work were obtained solely through the use of numerical and empirical methods, without the use complex analytical formulas that are typically used to describe similar 2D systems.This was necessary due to the additional analytical complexity required to evaluate and predict the maximum stresses in 3D systems and also since it allows one to reduce the relationships between the geometric parameters and stresses in 3D systems into relatively short formulas which are much more manageable and easier to use than analytical expressions.However, it must be emphasized that there are some limitations to this approach and that these empirical models may be used with maximum accuracy mainly to predict the stress concentration factor of structures which fall within the range of geometric parameters evaluated in this study, i.e.
within the range of t/d0 = 0.125 -1 and d0/d1 = 0.25 -1.This, of course, is a limitation of all empirical methods in general.For cases where the ligament length is extremely small (t/d0 → 0), the precision of these models decreases significantly as the steepness of the slope of the exponential equation reaches its maximum value and therefore, in order to accurately predict the stress concentration factor in these cases one must either conduct an additional numerical study of these cases (with separate mesh convergence testing criteria in order to cater for the extremely small ligament length) or a dual numerical-analytical approach using techniques such as asymptotic analysis [35].This is beyond the scope of this work and is an avenue which should be explored further in future studies on 3D skew hole systems.Moreover, it must be noted that these empirical models are also not designed to cater for skewed cross-bore systems, i.e. systems where the skewed holes intersect each other and thus t/d0 is negative.These cases must be considered as a completely separate entity and must be considered independently if the numerical-empirical approach is used to evaluate these systems.Nevertheless, despite these limitations, the approach used in this work is suitable to predict the maximum stresses found in a wide range of skew pressurized hole systems which include the main geometric cases commonly found in industrial settings.As shown in the plots presented previously, notwithstanding the relative reduction in complexity with respect to conventional analytical methods, the use of this approach does not result in a general lack of accuracy since there is good agreement between the empirical models and numerical results obtained.

Conclusion
In this work, we investigated the stress concentration factor found in pressurized non-aligned hole systems using Finite Element Analysis.A range of systems with various geometric parameters and different pressure-loading conditions were modelled and evaluated, and the numerical results were used to obtain empirical models which were used to create reference plots for the prediction of the stress concentration factors in such systems.The trends observed from the various changes in geometric parameters were also analysed using statistical and numerical methods; it was found that the t/d0 ratio is the most significant geometric parameter in determining the stress concentration factor.Furthermore, it was observed that in case of pressure applied in one hole only, the stress concentration factor is equal to the applied pressure when t/d0 ~ 1, indicating that the holes have negligible interactions between them and may be treated as a single hole in an infinitely large plane/block.The numerical results and the empirical models presented here also show good agreement with analytical results based on 2D models (see Appendix 1) and are expected to function as a useful tool for the pre-design and evaluation of stress concentration factors in non-aligned and non-concurrent pressurized hole systems.

Figure 1 :
Figure 1: A pictorial representation of the three 2D cases considered by Scire Mammano and Dragoni [16]: a) a pressurized hole inside an externally pressurized cylinder, b) a pressurized hole close to a pressurized boundary and c) two pressurized holes in close proximity within an infinite plane, where p1 and p0 represent two unequal pressures.Image d) shows a mock-up of a hydraulic manifold with non-aligned holes (image taken from[48]).

Figure 2 :
Figure 2: Images showing the simulated system from different orientations and indicating the geometric parameters used to define it.
The system was loaded by applying a uniform pressure, p, equivalent to 1 MPa on the internal surfaces of the cylindrical holes.Three load cases were set: Loadcase A both holes were subjected to an internal uniform pressure (p0 = p1 = 1 MPa), Loadcase B the first hole, with diameter d0, only was pressurized (p0 = 1 MPa, p1 = 0 MPa) and Loadcase C where the second diameter d1, only was pressurized (p0 = 0 MPa, p1 = 1 MPa).The simulations were solved linearly using the material properties of linear elastic structural steel, E = 210000 MPa and ν = 0.3.

Figure 3 :
Figure 3: Diagrams showing the constraints applied to the system and the cylindrical surfaces upon which a uniform pressure was applied.

Figure 4 :Figure 5 :
Figure 4: Plots showing the variation of the stress concentration factor with changes in θ for the various geometric systems studied here (t/d0 and d1/d0) under p0 = p1 = 1 MPa loading conditions.
f(x,y) is Kt, and x and y represent the geometric parameters θ (expressed in radians) and t/d0 respectively.The following empirical model was obtained, with an R 2 value of 0.994, where the constants are equal to: a = -0.290,b = 0.573, c = 0.050 and d = -0.843.

Figure 6 :
Figure 6: 2D plot showing the empirical model presented in Eq. 3 fitted over the data points.

5 exhibit
Kt values below 1.3.This value is very close to the unitary normalized stress value,

Figure 7 :
Figure 7: Plots showing the variation of the stress concentration factor with changes in θ for the various geometric systems studied here (t/d0 and d1/d0) using Loadcase B.
x,y,z) is Kt, and x, y and z represent the geometric parameters θ (expressed in radians) t/d0 and d1/d0 respectively and a, b, c, d, e and f are the constants that need to be found.The following empirical model was obtained:

Figure 8 :
Figure 8: Reference plots, including the FE data points, showing the curves obtained from the empirical model presented in Eq. 5.

FittingFigure 9 :
Figure 9: Plots showing the variation of the stress concentration factor with changes in θ for the various geometric systems studied here (t/d0 and d1/d0) using Loadcase C.The model presented in Eq. 6 was used for the curve-fitting of the results of this loadcase:

Figure 10 :
Figure 10: Reference plots, including the FE data points, showing the curves obtained from the empirical model presented in Eq. 7.
for the material considered, such as Tresca, Von Mises or Mohr criterion.