Formulization of Three-Dimensional Stress and Strain Fields at Elliptical Holes in Finite Thickness Plates

Stress raisers such as holes are inevitable in structures at which stress concentration occurs and the static as well as fatigue strength of the structures can be significantly weakened. Therefore, to accurately evaluate the stress concentration factor and stress fields at holes is of essential importance for structure design and service life prediction. Although stress and strain concentration and fields at holes in finite thickness plates strongly change with and along the thickness, manuals of stress concentration for engineering design are mainly based on two-dimensional theory and no explicit formula is available even for circular holes in finite thickness plates. Here we obtain for the first time a complete set of explicit formulae for stress and strain concentration factors and the out-of-plane constraint factor at circular as well as elliptical holes in finite thickness plates by integrating comprehensive three-dimensional finite element analyses and available theoretical solutions. The three-dimensional stress distributions ahead of holes can also be predicted by the obtained formulae. With their accuracy and the corresponding applicable range being analyzed and outlined in detail, the formulae can serve as an important fundamental solution for three-dimensional engineering structure design and guideline for developing three-dimensional analytical methods.


Introduction
Machines or structural components inevitably contain stress raisers such as holes that can serve as fatigue crack initiation points. Accurate evaluation of the stress concentration factor (SCF) and stress fields at the stress raisers becomes the prerequisite for structural design and service life prediction. Pioneer analytical solutions for SCFs at circular and elliptical holes in infinite thin plates were obtained by Kirsch [1] and Inglis [2]. Complex variable solutions for the same problem were obtained later [3][4][5]. Rational solutions were reported for holes and notches in infinite thin plates and empirical formulae were proposed for more practical cases [6][7][8][9][10][11][12]. However, these theoretical solutions and empirical formulae were mainly obtained within the two-dimensional (2D) frame.
In developing analytical three-dimensional (3D) solutions for SCFs at cavities and holes, outstanding efforts have been made since 1940s. Sadowsky and Sternberg [13,14]  near a general ellipsoidal cavity in an infinite body. Sternberg et al. [15] reported an approximate solution for the 3D stress distributions ahead of a circular hole in an infinite plate with arbitrary thickness, which was revisited with reference to uniaxial tensile loading [16]. A generalized plate theory considering the thickness effect was proposed by Kane and Mindlin [17], and has been widely applied recently [18][19][20][21] with the assumption of uniform through-the-thickness extensional strain in the thickness direction. Although the average thickness effect can be reflected in certain degree, the variations of the SCFs along the thickness direction cannot be predicted in the frame of generalized plate theory. An important achievement was made by Folias and Wang [22]. They presented a general analytical solution for circular holes in finite thickness plates on the basis of Navier's equation and managed to work out the 3D SCF along the thickness for the first time following an elaborate numerical process. They found that the SCF is sensitive to the plate thickness and Poisson's ratio and showed that the SCF attains its maximum on the mid-plane and the plane near the free surface for thin plate and thick plate, respectively. With the increasing capability of calculations, numerical analyses have been more and more widely used in investigating the SCFs and stress fields at holes and notches in finite thickness plates. Guo et al. [23][24][25][26][27] investigated the 3D stress fields in notched plates of finite thickness subjected to tensile loads. They found that the SCF along the 3D notch tip depends on notch configuration and the thickness to root radius ratio, and its value on the mid-plane can be 25% higher than that on the free surface and deviate about 8% from the 2D solution. Therefore, the use of 2D results cannot satisfy safe design of modern structures [28]. Later researchers [29][30][31] presented similar results and revealed that stress and strain concentration factors are not the same along the thickness even in linear elastic   7), the previous study [40] and the 2D FE analysis plates. The importance of 3D effect has been highlighted in some review articles [32,33]. Recently, Berto et al. [34] confirmed the existence of the stress field associated with the out-of-plane mode in U-notches, as well as circular and elliptical holes under in-plane shear loading. They also reported that SCFs are distinctly different between 2D and 3D models for elliptical holes under in-plane shear loading. Zappalorto and Carraro [35] published a useful engineering formula for the theoretical SCF of orthotropic notched plates under tension. This engineering formula needs, as inputs, the material elastic constants and the SCF of the corresponding isotropic case. These experimental and numerical evidences show that 3D stress and strain concentration and fields at holes and notches are of great significance for structural design. However, no explicit formula is available even for the distribution of SCF at a circular hole along the plate thickness in a finite thickness elastic isotropic plate yet. Explicit and practical expressions for 3D stress and strain fields are still of great significance in developing three-dimensional analytical methods and rapid evaluation for engineering practice.
In this article, we report for the first time a complete set of explicit formulae for 3D stress and strain concentration and fields at circular and elliptical holes in finite thickness plates subjected to uniaxial tensile load. The accuracy and applicable range of the formulae are outlined in detail to guide applications.

Definition
A large elastic isotropic plate with thickness B containing an elliptical hole, which is subjected to a far-end uniaxial tension, is considered here, as illustrated in Fig. 1a. The elliptical major axis is 2a, minor axis is 2b, and the elliptical shape factor is t = b/a for the hole. The corresponding curvature radius is ρ = b 2 /a = at 2 . The width of the plate is 2W . The height of the plate is 2H. The geometrical parameters of these specimens are as follows: a = 7 mm, H/a = W/a = 20, B/a = 0, 0.5, 1, 2, 3, 4, 5, 8, 10, 30 and t = 1, 0.8, 0.5, 0.3, 0.1. The origin of the Cartesian coordinate system (x, y, z) is at the center of the elliptical hole and the plate free surfaces are at z = ± B/2. For convenience, the corresponding cylindrical coordinate system (r, θ, z) is also adopted, as shown in Fig. 1a.
The SCF K t (r = a, θ = 0) is where σ net = σ y W/(W − a) is the mean stress of the net section on the ligament and σ y is the applied stress at the ends of the plate along the y direction. K t on the mid-plane of the plate (2z/B = 0) is denoted as (K t ) mp . Since the geometrical parameters 2H and 2W are sufficiently large compared with the elliptical major axis 2a, the SCF for plane stress state defined as (K t ) p−σ = σ yy /σ y can be calculated by the 2D finite element (FE) method. Three parameters are introduced to describe the 3D characteristics of the stress fields in the finite thickness plates. The first parameter is the in-plane stress ratio defined as The second parameter, the out-of-plane stress constraint factor, has been defined by Guo [36][37][38] as In an isotropic linear elastic cracked body, T z = 0 for the plane stress state and T z = ν for the plane strain state. For the finite thickness plate, 0 ≤ T z ≤ ν, where ν is Poisson's ratio of the material. The third parameter is the equivalent thickness B eq of the point on the crack front line empirically defined by She et al. [39]:

Finite Element (FE) Modeling
Due to symmetry in geometries, only one-eighth of each plate is modeled with finite elements (region x ≥ 0, y ≥ 0, 0 ≤ z ≤ B/2). All analyses are carried out using the software ABAQUS. As shown in Fig. 1c, the 3D mesh is constructed with the eight-node solid linear elements. In order to characterize the through-thickness stress gradient, forty element layers are divided along the half-thickness of the plate, and the thickness of the element layer is gradually increased from the free surface (2z/B = 1) to the mid-plane (2z/B = 0), with the bias ratio being 100 (t = 1, 0.8) and 1000 (t = 0.5, 0.3, 0.1) in mesh sizing control. Within each layer, the size of forty elements enclosing the surface of the hole gradually increases along the circumferential direction from θ = 90 • to θ = 0 • , where the bias ratio is 1000 in mesh sizing control. In the x-direction, the first forty elements starting from x = a are gradually enlarged with the bias ratio 1000 to reflect the in-plane stress gradient. In particular, the 2D plane model is adopted for B/a = 0. The 2D FE mesh is constructed with the four-node solid linear 2D elements, with the mesh sizing control strategy similar to the 3D model. These FE models represent a compromise between the required level of mesh refinement and the extensive computation time. The applied stress at the end of the plate σ y is 50 MPa. Poisson's ratio ν = 0.33 and Young's modulus E = 200 GPa.

SCF on the Mid-Plane
The SCF K t calculated by Eq. (1) is an important parameter scaling the extent of stress raisers around holes. The variations of the SCF on the mid-plane (2z/B = 0) (K t ) mp versus the dimensionless plate thickness B/a are plotted in Fig. 2a-e. (K t ) mp is not a monotonic function of dimensionless plate thickness (B/a) for circular hole (t = 1, a = ρ), which is similar to the conclusion by Li et al. [23]. However, its peak value (K t ) mp-max is at about B/a = 2, as shown in Fig. 2a, rather than at B/ρ = 3 ∼ 8. For t = 1 and 0.8, both (K t ) mp-max are at around B/a = 2. For t = 0.5 and 0.3, both (K t ) mp-max are at about B/a = 1. In addition, (K t ) mp-max is at about B/a = 0.5 for t = 0.1. With smaller t, the increasing rate becomes higher before reaching the peak value and the decreasing rate becomes higher after reaching the peak value. After decreasing from the peak value, (K t ) mp is stable at about B/a = 10. The difference between (K t ) mp-max and the value at B/a = 0 increases from 0.18 to 1.5 with the decrease of elliptical shape factor from t = 1 to t = 0.1. The SCF on the mid-plane at holes can be about 8% higher than the 2D solution. The relationships among (K t ) mp , B/a and t can be fitted by the FE results. For 0 ≤ B/a ≤ 30 and 0.1 ≤ t ≤ 1, (K t ) mp can be expressed as where q 2 (t) = 10.816t 3 − 6.3259t 2 − 15.3192t + 7.3886, q 3 (t) = 7.1503t 3 − 12.648t 2 + 4.3624t − 0.1671, q 4 (t) = −2.5714t 3 + 9.5321t 2 − 2.1436t + 0.5833.
(K t ) mp described with Eq. (5) is shown in Fig. 2a-e. As shown in Fig. 2f, the relative errors of Eq. (5) and the FE results are within ± 0.75% in the range of B/a = 1 ∼ 30 and t = 0.1 ∼ 1. When B/a approaches zero, q 1 (t) becomes the 2D SCF for infinite plates. As can be seen from Fig. 3, Eq. (7) is in good agreement with the 2D FE results as well as the formula proposed by Dai et al. [40], which is in the form of (K t ) p−σ = 1 + 2/t in the range of 0.01 ≤ t ≤ 1.  Figure 4 shows the variations of SCF K t normalized by (K t ) mp along the thickness from the midplane (2z/B = 0) to the free surface (2z/B = 1) in different elliptical holes for different B/a. The data points are the 3D FE results. For t = 1, 0.8, 0.5, 0.3, the minimum K t /(K t ) mp along the thickness is located in the free surface (2z/B = 1) for different B/a. However, the maximum K t /(K t ) mp along the thickness is not always located in the mid-plane (2z/B = 0) or the plane layer near the free surface. K t /(K t ) mp attains the maximum on the mid-plane for B/a no more than 3 at holes with t = 1 and 0.8. The maximum K t /(K t ) mp is on the mid-plane for B/a no more than 2 and 1 for t = 0.5 and 0.3, respectively. The value of B/a = 3, 2 or 1 is called the transition thickness (B * σ /a) proposed by Yang et al. [29]. They pointed out that the location of the maximum K t /(K t ) mp is not in the mid-plane or the plane layer near the free surface, but shifts gradually between the two, when the plate thickness exceeds the transition thickness. Figure 5 shows the variations of SCF K t along the thickness from the mid-plane to the free surface in different elliptical holes for different B/a. The complex distribution of K t reveals the inherent complexity of 3D problems. The data points are the 3D FE results, and the dash dot lines are the 2D FE results. The trend of K t changes slightly when K t is greater than the 2D solution. The K t near the free surface drops rapidly when it is smaller than the 2D solution. As shown in Fig. 5d, the K t near the mid-plane is not always larger than the 2D solution for t = 0.3. Therefore, the quantitative description of K t can reduce the risk of just using the SCF measured at the free surface or from the 2D results for engineering design. For different elliptical factors t, the relationship among the K t /(K t ) mp , 2z/B and B/a can be fitted by the FE results. For 1 ≤ B/a ≤ 10 and 0.3 ≤ t ≤ 1, the relationship can be expressed in the form of

SCF Along the Thickness
where The fitted value of K t /(K t ) mp described with Eq. (8) and the FE result of K t /(K t ) mp are plotted in Fig. 4. Equation (8) times Eq. (5) is Eq. (11) in the form of  11), it reads that in the interior of the plate the 3D SCF at the circular hole with t = 1 is within 6% above the 2D solution, while the 3D SCF at the circular hole with t = 1 can be 16.5% lower than the 2D solution at the free surface. For the elliptical hole with t = 0.3, the deviations of the 3D results from the 2D solutions are within 5.4% in the interior but can be 29.5% at the free surface. Figure 6 shows the relative error of the fitted value of K t and the FE result of K t versus 2z/B for different B/a changing from 1 to 10. The relative error is within ± 1.25% in the range of 0 ≤ 2z/B ≤ 0.95 for elliptical factor t = 1; meanwhile, it is within ± 1.25% in the range of 0 ≤ 2z/B ≤ 0.86 and 0 ≤ 2z/B ≤ 0.96 for t = 0.8 and 0.3, respectively. When it comes to t = 1 and 0.3, the fitted value of K t has the accuracy of over 5% just near the free surface only for B/a = 10 and 8. In particular, the relative error is within ± 1% in the range of 0 ≤ 2z/B ≤ 1 for t = 0.5. Figure 7 shows 3D distributions of the fitted value of K t versus B/a and 2z/B in different elliptical holes. The surface represents Eq. (11), and the solid lines are plotted according to the FE results. In particular, the FE results at 2z/B = 1 for B/a ranging from 0.5 to 10 are marked as black dots at the ends of the solid lines. Accuracy of Eq. (11) decreases near the free surface for t = 1, 0.8, 0.5 when B/a = 0.5. Equation (11) is in good agreement with the FE results for t = 0.3.

In-Plane Stress Ratio T x and In-Plane Stress Components σ y y , σ xx at Elliptical Holes
Previous studies on constructing approximate formulae for stress distribution ahead of a notch tip have been intensively investigated [8,41,42]. The best-known expression is the classical solution [5] to the stress fields in an infinite plate containing an elliptical hole under uniaxial loading condition. This solution is in the form of Eqs. (12)-(14), when SCF is introduced in 2D cases: It is easy to derive from Eqs. (12) and (13) that It should be pointed out that Eq. (14) is valid in a large range of (r −a)/ρ for the circular hole, but only within a small distance near the elliptical hole or notch root such as (r − a)/ρ ≤ 0.3 [23]. Analogous to the forms of Eqs. (12)- (14), Eqs. (15)- (17) are obtained for elliptical holes with t ranging from 1 to where (14) when t = 1. Figure 8 shows the comparison between formulae and the 2D FE results of in-plane stress components and in-plane stress ratio versus (r − a)/(at 2 ). Equations (15) and (16) are in good agreement with the 2D FE results in a large range of 0 ≤ (r − a)/(at 2 ) ≤ 100. As can be seen from parts (a) and (b) of Fig. 8, parts of the valid scope of (r − a)/(at 2 ) are chosen to distinguish the distributions of FE results as well as the fitted results at different elliptical holes (t = 1, 0.8, 0.5, 0.3). The in-plane stress ratio T x in Fig. 8c increases from 0 to its peak value and then decreases gradually to 0 with increasing (r − a)/(at 2 ). The smaller is the elliptical factor t, the larger is the peak value of in-plane stress ratio T x . There is an overlap between the 2D FE results for different t in the range of 0 ≤ (r −a)/(at 2 ) ≤ 0.3, which can be described by Eq. (14) for different notch configurations. As shown in Fig. 8c, the fitted results described with Eq. (17) are consistent with the FE results. Based on Eq. (16) in 2D cases, we introduce the 3D SCF to describe the in-plane stress component σ yy | θ=0 /σ y on different plane layers in finite thickness plates, which can be expressed as where K t (2z/B, B/a, t) has been derived in Sect. 3.1.2 and (K t ) p−σ = σ yy /σ y is in Sect. 2.1.
To distinguish the distributions of σ yy | θ=0 /σ y on different plane layers, parts of the valid scope of (r − a)/(at 2 ) are chosen in Figs. 12, 13 and 14. Figures 12, 13 and 14 show the distributions of σ yy | θ=0 on different plane layers for dimensionless thickness B/a = 1, 5, 10, respectively. Equation (18) can predict the distributions of σ yy | θ=0 very well on different plane layers except for the free surface (2z/B = 1) at t = 1, 0.3 for B/a = 10. The prediction with low accuracy attributes to the relatively large error of Eq. (11) in some cases including but not limited to the description of K t at the free surface at t = 1, 0.3 for B/a = 10.
Equation (18)   , while the dash lines represent the 2D formula of Eq. (15). In particular, the distributions of σ xx | θ=0 /σ y can be well predicted with Eq. (19) even at the free surface.

T z and σ zz at the Root of Holes
The local 3D stress field becomes very strong when an elliptical hole degenerates from a circular hole (t = 1) to a crack (t = 0). The out-of-plane stress constraint factor T z defined by Eq. (3) has been extensively studied to characterize the out-of-plane stress component near the crack front. The saturated value of T z /ν is reported as 2/3 for the circular hole [14]. However, no explicit formula of saturated value of T z /ν versus plate thickness for elliptical holes has been reported. By fitting of the 3D FE analysis results with sufficiently high accuracy, the out-of-plane stress constraint factor normalized by Poisson's ratio at the root of holes, T z0 /ν, changing with and along the finite thickness at elliptical holes can be described as where a(t) = − 2.2124t 3 + 4.5951t 2 − 3.752t + 2.0118, b(t) = − 2.386t 3 + 5.2084t 2 − 4.6204t + 2.1116 and c(t) = − 1.0587t 3 + 6.9335t 2 + 0.1583t + 0.0333. Figure 18 shows the variations of the out-of-plane stress constraint factor normalized by Poisson's ratio at the root of holes of the mid-plane, T z0 /ν, versus B/a for elliptical holes. T z0 /ν increases with the dimensionless plate thickness B/a and tends to be stable when B/a exceeds a certain value, which is usually expressed as the sufficiently large thickness. T z0 /ν becomes stable when B/a attains around 10 which is almost independent of the elliptical shape factor t. Equation (20) is in good agreement with the FE results on the mid-plane (2z/B = 0). Figure 19 shows T z0 /ν on the mid-plane (2z/B = 0) versus the elliptical factor t in a finite thickness plate with B/a = 10. The smaller is the elliptical shape factor t, the larger is the saturated value of T z0 /ν. This is supported by the FE results in the literature [23,30].
The variations of the out-of-plane stress constraint factor normalized by Poisson's ratio at the root of holes, T z0 /ν, versus 2z/B for different elliptical shape factors t are plotted in Fig. 20. T z0 /ν attains its maximum on the mid-plane and decreases monotonically from the mid-plane (2z/B = 0) to the free surface (2z/B = 1).
With respect to B eq /a rather than 2z/B, T z0 /ν can be well predicted with one curve for a given elliptical factor t shown in Fig. 21. Figure 22 shows the differences between the formula values and the FE values of T z0 normalized by Poisson's ratio versus B eq /a for different elliptical shape factors t. Most of the error values are in the range of ± 2.5%, and the maximum error is smaller than 5%. According to the definition of T z and T x , the out-of-plane stress component can be expressed as When r = a and θ = 0, Eq. (23) can be obtained as According to Eqs. (11) and (20), Eq. (23) can be expressed as Figure 23 shows the variations of the out-of-plane stress component σ zz normalized by the applied stress at the ends of the plate σ y versus 2z/B for different elliptical shape factors t. The out-of-plane stress component σ zz attains its maximum on the mid-plane and decreases monotonically from the mid-plane (2z/B = 0) to the free surface (2z/B = 1). The FE results can be well predicted with Eq. (24). Figure 24 shows the differences between the fitted values and the FE results of σ zz normalized by the applied stress at the ends of the plate σ y versus 2z/B for different elliptical shape factors t. Most of the error values are in the range of ± 5% for B/a ranging from 1 to 10 for t = 1, 0.8 and 0.5. For t = 0.3, the error value is between minus 4% and minus 10% with 2z/B ranging from 0 to 0.6, and is within ± 5% with 2z/B ranging from 0.6 to 1.

T z and σ zz in 3D Affected Zone
As is well known, a strong 3D stress field is observed within the radius of half plate thickness from the crack tip in a sufficiently thin plate. For U-notch, the 3D affected zone (T z > 0) on the mid-plane where (T z ) mp is T z on the mid-plane and (T z0 ) mp is T z at the root of holes of the mid-plane. Equation (25) needs verifying for larger B/a and different elliptical holes. By fitting of the 3D FE results for B/a = 10, when the concept of equivalent thickness B eq is introduced, we approximate the distribution of T z within the 3D affected zone on different plane layers for elliptical holes by where α(t) = −0.315t 2 + 1.8211t + 2.4165, and β(t) = 1.4745t 2 − 3.4915t + 2.8012. m(t) and n(t) have been proposed in Eq. (16). The maximum root-mean-square error is 0.02712 for different elliptical holes (t = 1, 0.8, 0.5, 0.3). Figure 25 shows the variations of the out-of-plane constraint factor T z normalized by T z0 , the out-ofplane constraint factor at the root of holes, versus (r − a)/B on different plane layers for dimensionless thickness B/a = 10. T z /T z0 is maximal at the root of holes and decays with the distance from the root of holes. Equations (25) and (26) are plotted in Fig. 25 when 2z/B = 0. Equation (25) can well predict the 3D FE results only on the mid-plane at the circular hole (t = 1) for B/a = 10. Figure 26 shows the variations of T z normalized by T z0 versus (r − a)/B eq on different plane layers for dimensionless thickness B/a = 10. The 3D FE results can be mostly predicted with one curve described as Eq.  Fig. 28, the 3D affected zone is approximately limited to a distance of 0.375B eq , 0.4B eq , 0.45B eq , 0.5B eq from the root of holes for t = 1, 0.8, 0.5, 0.3, respectively. T z /T z0 in the 3D affected zone can be mostly predicted with Eq. (25), in which B is replaced with B eq on different plane layers except for the one near the free surface. As shown in Fig. 29, the 3D affected zone is approximately limited to a distance of 0.5B eq from the root of holes for B/a = 5. The size of 3D affected zone is almost independent of the elliptical shape factor t ranging from 1 to 0.3.

Stress Triaxiality R σ
The stress triaxiality which dominates ductile fracture is defined as where σ m is hydrostatic stress and σ e is the von Mises equivalent stress Neglecting σ xy , σ yz and σ zx in the 3D affected zone, the stress triaxiality R σ for θ = 0 can be expressed as Substituting Eqs. (17), (20), and (26) into Eq. (32) yields, Figure 34 shows the variations of R σ at the root of holes (r = a) versus 2z/B for different B/a. R σ decreases monotonically from the mid-plane to the free surface along the thickness direction for given plate dimensionless thickness B/a and elliptical factor t. For given plate dimensionless thickness B/a and plane layer 2z/B, the smaller is the elliptical factor t, the larger is R σ . Equation (32) can predict the variations of R σ at the root of holes for different B/a in different elliptical holes (t = 1, 0.8, 0.5, 0.3) with the maximum relative error smaller than 4%. Figures 35 and 36 show the distributions of R σ on different plane layers for dimensionless thickness B/a = 1 and 10, respectively. With the distance to the root of holes increasing, R σ increases to its peak value and then decreases for different plane layers in different thickness plates containing an elliptical hole. With smaller elliptical shape factor t, the location of maximum R σ is much closer to the root of holes. As shown in Fig. 35, the distributions of R σ on different plane layers are very similar for t = 1 and 0.8. Equation (33) can predict the 3D FE results with adequate accuracy for engineering design at elliptical holes (t = 1, 0.8, 0.5). The predicted value of R σ is greater than the FE results for t = 0.5 and 0.3. R σ attains its peak value when the distance to the root of the hole increases to around 0.52a, 0.5a, 0.15a, 0.12a for t = 1, 0.8, 0.5, 0.3, respectively. As can be seen from Fig. 36, R σ attains its peak value when the distance to the root of holes increases to around a, 0.75a, 0.5a, 0.25a for t = 1, 0.8, 0.5, 0.3, respectively. Equation (34) cannot provide accurate prediction of R σ in quantitative terms. Even so, Eq. (34) can describe the general trend of increasing to a peak value and then decreasing. The location of maximum R σ of Eq. (34) is closer to the root of the hole than the FE results. The distance between the predicted location and the actual location is around 0.25a.

Strain Concentration Factor Along the Thickness
The strain concentration factor K ε (r = a, θ = 0) is where ε net = σ net /E is the mean strain of the net section on the ligament. For the 3D stress state of elastic body, the strain components can be denoted as ε yy = 1 E [σ yy − ν(σ zz + σ xx )]. According to Eqs. (1)-(3), the following expression can be obtained The in-plane stress ratio T x is zero along the root of holes. From Eq. (36), the distribution of K ε /K t along the root of holes can be expressed as where T z0 is described in Eq. (20). The strain concentration factor is distinctly different from the stress concentration factor along the thickness of the root of holes in finite thickness plates except at the free surface, as shown in Fig. 37. The K ε /K t at the mid-plane is smaller for larger thickness and smaller elliptical shape factor t. Similar distributions of K ε /K t along the root of the holes for different B/a as well as t were reported [30]. However, no explicit formula is available to describe the relationship between the strain concentration factor and the stress concentration factor yet. The distributions of strain concentration factor over stress concentration factor K ε /K t can be described with Eq. (37). The maximum relative error between Eq. (37) and the FE results is smaller than 0.8%.

Strain Energy Density Along the Thickness
In the case of linear elastic bodies, the elastic strain energy densities for uniaxial and multiaxial stress states can be expressed as Eqs. (38) and (39), respectively.
where R σ are shown in Eqs. (32)- (34). The variations of U/U 0 at the root of holes versus 2z/B for different elliptical shape factors t are shown in Fig. 38. U/U 0 decreases monotonically from the mid-plane to the free surface along the thickness direction for given plate dimensionless thickness B/a and elliptical factor t. For given plate dimensionless thickness B/a and plane layer 2z/B, the smaller is the elliptical factor t, the larger is U/U 0 . Equation (40) can predict the variations of U/U 0 at the root of holes for different B/a in different elliptical holes (t = 1, 0.8, 0.5, 0.3) with the maximum relative error smaller than 2%.  The stress and strain at the root of holes are necessities for predicting the fatigue life with traditional local strain approaches. An engineering method for strain hardening materials has been developed to predict the elastic-plastic stress-strain distribution ahead of the root of holes based on the elastic solutions [43]. Neuber's rule and the equivalent strain energy density (ESED) method were also adopted in this engineering method. The explicit formulae for stress and strain concentration factors as well as stress fields mentioned above can be used to modify Neuber's rule and the ESED method.

Conclusions
Based on comprehensive three-dimensional finite element analyses and the existing theoretical research, the three-dimensional distributions of stress and strain concentration factors at circular as well as elliptical holes and the stress fields near the holes are fitted into a set of explicit empirical formulae in finite thickness plates. Necessary error analysis and scope of validity are conducted for these explicit formulae. Equations (5) and (7) are also in good agreement with the available results in the literature. Equation (5) shows that SCF on the mid-plane at holes with different elliptical shape factors t = 1 ∼ 0.1 can be about 8% higher than the corresponding 2D solution. Equation (7) shows that SCF (K t ) mp on the mid-plane (2z/B = 0) at elliptical holes becomes the 2D SCF for infinite plates when B/a approaches zero. The 3D SCF K t changes above the 2D result with deviation within 6% in the interior of the plate, while it drops rapidly near the free surface and can be 29.5% below the 2D result for t = 0.3 to 1 and B/a = 1 to 10. The three-dimensional distributions of in-plane as well as out-of-plane stresses ahead of the holes can also be predicted by the presented formulae efficiently. Therefore, the obtained formulae can serve as useful basic explicit solutions for engineering structure design and guideline for developing three-dimensional analytical methods. The following conclusions can be helpful for application of the formulae.   (3) In the range of 0 ≤ B/a ≤ 30, the out-of-plane stress constraint factor normalized by Poisson's ratio, T z /ν, can be described by Eq. (18), changing from the mid-plane to the free surface at the root of elliptical holes (0.1 ≤ t ≤ 1). With respect to B eq /a rather than 2z/B, T z /ν can be well predicted with one curve for a given elliptical factor t. Most of the error values are in the range of ± 2.5% and the maximum error is within 5%. In the range of 1 ≤ B/a ≤ 10, Eq. (22) is obtained to characterize the out-of-plane stress component σ zz at the root of holes normalized by the applied stress at the ends of the plate σ y with errors in the range of ± 5% for elliptical holes (t = 1, 0.8, 0.5) and − 10% ∼ + 5% for t = 0.3. (4) With increasing distance to the root of the hole, the stress triaxiality R σ increases to its peak value and then decreases for different plane layers in different thickness plates containing an elliptical hole. Equation (32) can predict the variations of the stress triaxiality R σ at the root of elliptical hole for different B/a in different elliptical holes (t = 1, 0.8, 0.5, 0.3) with the error within 4%.

Stress Fields Within the 3D Affected Zone
(1) A set of formulae, Eqs. (15)- (19), is obtained for holes with elliptical shape factor t ranging from 1 to 0.3. The 2D solution of Eq. (17) can describe the in-plane stress ratio T x on different plane layers in finite thickness plates. The in-plane stress components σ xx and σ yy can be predicted very well by Eqs. (18) and (19) on different plane layers. The 3D affected zone can be described as the multiple of B eq on different plane layers. (2) The distributions of the out-of-plane constraint factor T z normalized by T z0 , the out-of-plane constraint factor at the root of holes, can be described with Eqs. (25) and (26) for B/a = 1 and 10, respectively. The out-of-plane stress component σ zz can be well predicted by Eqs. (28) and (27) for B/a = 1 and 10, respectively. (3) When B/a = 1, R σ attains its peak value when the distance to the root of the hole increases to around 0.52a, 0.5a, 0.15a, 0.12a for t = 1, 0.8, 0.5, 0.3, respectively. Equation (33) can predict the 3D finite element results of R σ with adequate accuracy for engineering design at elliptical holes (t = 1, 0.8, 0.5). When B/a = 10, R σ attains its peak value when the distance to the root of the hole increases to around a, 0.75a, 0.5a, 0.25a for t = 1, 0.8, 0.5, 0.3, respectively. Equation (34) can describe the general trend of increasing to the peak value of stress triaxiality R σ and then decreasing. However, the location of maximum R σ by Eq. (34) is closer to the root of the hole than the finite element results. The distance between the predicted location and the actual location is around 0.25a. See Table 1.