Three-Dimensional Stress Fields in Thick Orthotropic Plates with Sharply Curved Notches under In-Plane and Out-of-Plane Shear

In this paper, an analytical solution for the stress fields in the close neighbourhoods of radiused notches in thick orthotropic plates under shear loading and twisting is provided. In the first step, the equations of the three-dimensional theory of elasticity are successfully reduced to two uncoupled equations in two-dimensional space. Later, the 3D stress field solution for orthotropic plates with radiused notches is presented and its degree of accuracy is discussed by comparing theoretical results and numerical data from 3D FE analyses. The solution proposed can be satisfactorily used to characterise the stress field in plates made with polymeric composite materials, such as fibre-reinforced polymers and natural composites.


Introduction
The regions close to geometrical variations, i.e., notches, holes and cutouts, are characterized by a significant concentration of the stress fields, which are responsible for crack formation under static and cyclic loading conditions.
Predicting the fatigue crack nucleation or the static failure of notched solids is of primary concern for mechanical designers, and the formulation of effective strength criteria stems from the accurate knowledge of the local stress fields, thus justifying the significant attention paid by different researchers on this topic over the years (see, among the others, Refs. [1][2][3][4][5][6][7][8][9][10][11][12] and references quoted therein).
All the above-mentioned papers refer to two-dimensional analyses, even if the mechanical components are three-dimensional by nature, and the number of analytical works in the literature devoted to 3D stress fields is very limited, especially in the case of orthotropic materials.
Dealing with this topic, it is worth mentioning the scientific results by Prabhu and Lambros published in 2000 [13], where a numerical investigation on the 3D stress distribution in a cracked orthotropic plate is carried out, and by Cheng et al. [14] who, in 2015, studied the V-notch problem in an orthotropic solid under a generalised plane deformation state. Differently, Pageau and Biggers in 1996 [15] developed a 3D finite element formulation to derive the singular stress state near material and geometric discontinuities in orthotropic and anisotropic media, whereas Choi and Folias in 1994 [16] analysed the three-dimensional stress fields in a composite plate weakened by a hole.
On parallel tracks, in the early 2000s, Kotousov and Wang [17,18] studied transversely isotropic composite plates with holes or inclusions, proposing a generalised plane-strain theory, whilst later on, Zappalorto and Carraro studied thick anisotropic plates with pointed V-notches [19,20], and in 2017, developed a new three-dimensional theory for anisotropic blunt cracks [21].
Moving to recent progresses in the field, worthy of mention are the results of Shi et al. [22], who developed a method for designing thin cylindrical shells for aerospace applications

•
To provide evidence that the 3D solution derived by Zappalorto and Carraro [19] for pointed notches can be extended also to orthotropic plates with holes or lateral radiused notches with any notch opening angle, under the hypothesis of a sufficiently small notch tip radius; • To show that, on the basis of the plane solution, stress components σ xx , σ yy and τ xy in a thick 3D anisotropic notched plate (i.e., the in-plane stress fields) can be accurately determined, whereas out-of-plane shear stresses, τ xz and τ yz , can be assessed using the pure antiplane shear solution; • To also document the presence of coupled modes for orthotropic thick plates weakened by holes or lateral notches with any notch opening angle, as reported in other research articles concerning isotropic components.
Finally, it is worth mentioning that the obtained solution is very useful to analyse the stress fields in polymeric composites, such as carbon-fibre-reinforced epoxy laminates, short fibre-reinforced polymers or natural composites.

Basic Field Equation
Consider a thick orthotropic plate weakened by a radiused notch. With reference to a Cartesian orthogonal system centred on the symmetry plane of the plate (see Figure 1 or Figure 2), the stress-strain relationships can be written as: where S ij are conventional compliance coefficients or, equivalently: where C ij are stiffness coefficients.     Figure 1. Elliptic hole in a three-dimensional plate under in-plane shear (a) or twisting (b). In (c) a schematic of elliptic hole is reported together with the dimensions used in the mathematical treatise and in the results' discussion.     Figure 2. Rounded V-shaped (hyperbolic) notches in a three-dimensional plate under in-plane shear (a) or twisting (b). In (c) a schematic of elliptic hole is reported together with the dimensions used in the mathematical treatise and in the results' discussion. Suppose now that displacement fields near the notch tip could be written in the following form by means of separated variables [19,28]: Invoking strain definitions, one obtains: Under the hypothesis that the tip radius of the radiused notch is small enough to assure that f (z) × u(x, y), f (z) × v(x, y) and g (z) × w(x, y) can be disregarded when compared to g(z) · ∂w(x, y)/∂x and g(z) · ∂w(x, y)/∂y, the strain field can be written as [19][20][21]: As a result of Equation (5), a state of "quasi" plane strain is present in the stress concentration areas near the notch tip; indeed, the through-the-thickness normal strain can be disregarded, whereas the out-of-plane shear strains cannot. Equation (5) further gives: ∂u(x,y) ∂x Moreover, Equation (5) combined with Equation (1) gives the possibility to formulate the stress components according to the following expressions: Consider now the equilibrium equation in the z direction: Substituting Equation (6) into Equation (8) results in: and since f (z) × u(x, y) and f (z) × v(x, y) can be disregarded if compared to g(z) · ∂w(x, y)/∂x and g(z) · ∂w(x, y)/∂y, Equation (9) further simplifies to give: where C ij are stiffness coefficients, or, more explicitly: Concerning the equilibrium equations in the x-direction and in the y-direction, the Airy stress function, φ, defined in the (x,y) plane, guarantees that such equilibrium is preserved [19][20][21]: Substituting Equation (12) into Equation (1), further invoking the compatibility equation in the (x,y) plane gives the possibility to formulate the following fourth-order differential equation: where: Equations (10) and (13) must be both valid at the same time. Accordingly, by solving the following two differential equations (which are uncoupled) it is possible to simplify the equations that describe the three-dimensional elastic fields for a plate with a notch: One should note that Equation (15) is linked to the plane strain problem of the orthotropic theory of elasticity, whereas Equation (16) is the equation defining the solution of the antiplane one.
The practical consequence of Equations (15) and Equation (16) is that the application of an out-of-plane shear stress, expressed by a non-zero w function, induces a local non-zero Airy φ function due to three-dimensional effects. At the same time, the application of any kind of in-plane loading (such as bending, tension or in-plane shear) generates a local non-zero w function. This coupled behaviour is not accounted for by the standard plane strain or plain stress solutions.
In particular, load conditions applied to a thick plate resulting in a skew-symmetric φ function (mode 2) create, due to three-dimensional effects, an antisymmetric w function (mode 3) in the close neighbourhood of the notch tip.
In Equations (17) and (18), µ 1 and µ 2 are unequal imaginary numbers defined as: and represent the roots of the following characteristic equation [7,10,30]: under the condition that (2B 12 + B 66 ) 2 ≥ 4B 11 B 22 . Differently, z 1 and z 2 are complex variables which are defined as: Equation (18) can be conveniently re-written as [10]: where the auxiliary angular functions introduced, k ij , m ij and n ij , read as follows:

Solution for the Quasi-Harmonic Equation (Out-of-Plane Shear Stresses)
Out-of-plane shear stress fields in the considered thick orthotropic plate can be written in terms of a complex function [31]: in cartesian coordinates, or: in polar coordinates, where ϕ 3 is a proper complex function to be chosen depending on the specific notch geometry under consideration, z 3 = x + µ 3 y, and:

Elliptical Hole in a Thick Plate under Shear
Consider a thick plate with an elliptical hole under shear or torsion (see Figure 1). The in-plane stress fields can be determined according to the solution proposed by Savin [32]. Re-arranging Savin equations allows the explicit solution to this problem to be obtained, according to which the stress field is: where τ Max xy (z) is, for a given z value, the maximum shear stress occurring along the notch bisector line at a certain distance, x Max -a, from the notch tip, whereas: Along the notch bisector line, Equation (32) simplifies to give: where: The out-of-plane stress fields, instead, can be derived starting from the solution obtained by Zappalorto and Salviato [31] for the pure antiplane shear loadings, which can be re-written as: where: and τ Max zy (z) is, for a given z value, the maximum shear stress occurring at the notch tip. Along the notch bisector line, Equation (39) simplifies to give:

Lateral Radiused Notch under Shear
Consider a thick plate with two lateral radiused notches with a generic notch opening angle (2α) and subjected to shear or torsion (see Figure 2). According to the previous literature (see, among the others, [33], and references reported therein), the notch boundary can be described with a hyperbola-like curve with parametric equations: where: and ρ is the notch root radius. The in-plane stress fields can be determined according to the solution proposed by Pastrello et al. [34]: where, denoting with x and y the distances from the notch tip: and λ 2 is a linear elastic eigenvalue to be determined by solving the following equation: where γ = π − α. Parameters t 2 , µ 2 , ζ 2 , χ 12 , χ 21 , χ 22 and χ 23 depend on the notch geometry and the material properties and can be determined according to the procedure proposed in Ref. [34]. Moreover, one should note that parameter A(z) in Equations (45)-(47) can be expressed, for a given z value, as a function of the maximum shear stress along the notch bisector, the maximum principal stress along the notched edge or as a function of a mode 2 generalised stress intensity factor for that given plane (z value).
The out-of-plane shear stress components can be determined according to the solution proposed in Ref. [31], according to which the stress field can be written as follows: where τ MAx zθ (z) is, for a given z value, the maximum shear stress occurring at the notch tip, whereas: In Equation (51), x and y are the distances from the notch tip in the x and y directions, respectively.
Moreover, the following equations hold valid: and Along the notch bisector line Equation (50) simplifies to give: One should note that, in the case of an isotropic material, β 3 = 1 and, accordingly, λ 3 = π 2γ , according to [5,35].

Discussion and Results
In order to support the theoretical results previously described, a set of numerical simulations has been carried out. In more details, three-dimensional elastic FE analyses were performed using three different material systems, of which the elastic properties are summarised in Table 1. They represent typical properties of some polymeric composites: • Material 1 represents a unidirectional carbon-fibre-reinforced epoxy laminate with the fibres oriented in the direction of the notch bisector; • Material 2 represents the same material with fibres oriented in the direction normal to the notch bisector; • Material 3 represents a quasi-isotropic carbon-fibre-reinforced epoxy laminate (e.g., [(0/±45/90) n ] S ). Numerical investigations were carried out on thick notched discs (see Figure 3) by applying on their external surface a pure mode 2 or a pure mode 3 displacement field (see the expressions in Appendix A), in order to generate an in-plane or out-of-plane stress field in the notched solids. This modelling strategy was already used by Berto et al. [25] in order to ease the analysis of the local behaviour of a generic 3D solid. summarised in Table 1. They represent typical properties of some polymeric composites:


Material 1 represents a unidirectional carbon-fibre-reinforced epoxy laminate with the fibres oriented in the direction of the notch bisector;  Material 2 represents the same material with fibres oriented in the direction normal to the notch bisector;  Material 3 represents a quasi-isotropic carbon-fibre-reinforced epoxy laminate (e.g., [(0/±45/90)n]S). Numerical investigations were carried out on thick notched discs (see Figure 3) by applying on their external surface a pure mode 2 or a pure mode 3 displacement field (see the expressions in Appendix A), in order to generate an in-plane or out-of-plane stress field in the notched solids. This modelling strategy was already used by Berto et al. [25] in order to ease the analysis of the local behaviour of a generic 3D solid. The following notch geometries were considered:


Semi-elliptical notches with notch depth a = 5 mm and different notch root radii (  = 0.001 mm, 0.01 mm, 0.1 mm and 1 mm); in this case, the radius of the disc was 5 mm and various thicknesses were considered (t = 1 mm, 5 mm and 10 mm);  Hyperbolic notches with a notch depth a = 10 mm, a notch root radius  = 0.01 mm and various notch opening angles (see Tables 2 and 3 for the stress field parameters associated with the cases analysed).  The following notch geometries were considered: • Semi-elliptical notches with notch depth a = 5 mm and different notch root radii (ρ = 0.001 mm, 0.01 mm, 0.1 mm and 1 mm); in this case, the radius of the disc was 5 mm and various thicknesses were considered (t = 1 mm, 5 mm and 10 mm); • Hyperbolic notches with a notch depth a = 10 mm, a notch root radius ρ = 0.01 mm and various notch opening angles (see Tables 2 and 3 for the stress field parameters associated with the cases analysed).  FE analyses were carried out with the commercial FE code ANSYS, release 19.5, using 20 nodes SOLID186 brick elements with reduced integration and pure displacement formulation.
The loads were applied as nodal displacements on the circular boundary of the disc (see Figure 3a). In more detail, the nodal coordinates of the boundaries were collected in ANSYS APDL, and, by means of an external script, the displacements to be applied were evaluated according to the expressions reported in Appendix A. Eventually, another APDL script was used to apply the nodal displacements to the disc boundary. For simulating mode 2, the symmetry of the model was exploited in order to use only one-quarter of the disc. Symmetry boundary conditions were applied to the face visible in Figure 3a, at a z-coordinate equal to half the disc thickness (i.e., the disc midplane). Anti-symmetry was applied to the nodes lying on the other symmetry plane. For simulating mode 3, only the anti-symmetry boundary conditions were used. Hence, the model needed was half of a full disc.
An important aspect to be accounted for in the definition of the FE model parameters is the elastic property of the plate material. In fact, depending on the material orientation, the mesh requires a different number of elements to reach a convergent solution. To minimize the computational cost of the model, the number of elements and the spacing ratio, both in the in-plane and through-the-thickness directions, required a multi-variable optimization. For each model, the parameters were modified until the maximum value of τ Max zy (z) at notch tip showed a variation of less than 1%. The choice of that output variable for the convergence analysis was made under the rationale that the most critical zone for mode 2 to mode 3 coupling is near the notch tip.
Further details about the numerical validation (e.g., mesh construction for each geometry) can be found in the Supplementary Material retrievable online.
Initially, the attention was focused on discs with a thickness t = 1 mm loaded under pure mode 2 and weakened by elliptical notches with a root radius ρ= 0.001 mm. Results related to the plane z/t = 0.4 (where z is the distance from the mid-plane) are presented in Figures 4-8. In particular, in Figure 4, the normal in-plane stress component tangent to the notch edge, σ vv , is plotted against the polar angle θ, measured with respect to the ellipse foci. Differently, the shear stress τ xy , as evaluated along the notch bisector line, is presented in Figure 5 as a function of the distance from the notch tip. As evident, in both cases, the results from the three-dimensional numerical analyses perfectly agree with the two-dimensional solution, Equations (32) and (37).
In Figures 6 and 7, the attention is instead focused on out-of-plane shear stresses induced by three-dimensional effects. As evident, these stress components can be accurately predicted using the solutions derived for the antiplane deformation, Equations (39) and (42). Figure 8 makes it evident that, considering several z/t values, the induced out-of-plane shear stresses remain self-similar for various z/t values, except when very close to the free surface of the disc. Accordingly, the value of z/t = 0.4 can be regarded as representative of any z value at a sufficient distance from the free edge of the disc. The particular choice of z/t = 0.4 was made to obtain comparable magnitudes between the induced stresses and the main stresses generated by the imposed loading conditions in the numerical simulations.           In Figures 6 and 7, the attention is instead focused on out-of-plane shear stresses induced by three-dimensional effects. As evident, these stress components can be accurately predicted using the solutions derived for the antiplane deformation, Equation (39) and Equation (42). Figure 8 makes it evident that, considering several z/t values, the induced out-of-plane shear stresses remain self-similar for various z/t values, except when very close to the free surface of the disc. Accordingly, the value of z/t = 0.4 can be regarded as representative of any z value at a sufficient distance from the free edge of the disc. The particular choice of z/t = 0.4 was made to obtain comparable magnitudes between the induced stresses and the main stresses generated by the imposed loading conditions in the numerical simulations.
Eventually, considering different values for the disc thickness (t = 1 mm, 5 mm and Eventually, considering different values for the disc thickness (t = 1 mm, 5 mm and 10 mm) the trend of the maximum in-plane shear stress, τ Max xy (z) and of the maximum out-of-plane shear stress τ Max zy (z) are reported in Figure 9a,b, respectively, where it can be noted that:

•
The maximum in-plane shear stress τ Max xy (z) remains constant for most of the plane thickness, and significantly increases as approaching the free surface of the plate (Figure 9a); • The maximum out-of-plane shear stress τ Max zy (z) has an almost linear trend up to z/t around 0.3. When approaching the free surface of the disc the trend becomes strongly nonlinear; • Increasing the disc thickness, the intensity of the induced out-of-plane shear stresses significantly increases (Figure 9b), whereas the influence of the thickness on τ Max xy (z) is limited (Figure 9a). Figures 10 and 11 present the effects of different notch root radii in terms of out-ofplane shear stresses. In particular, it can be noted from Figure 10 that, when increasing the notch radius:

•
The accuracy of Equation (42) decreases. Recalling the analytical treatise presented in Section 2, this is due to the terms f (z) × u(x, y) and f (z) × v(x, y), which are more and more relevant when decreasing the value of the notch root radius, when compared to g(z) · ∂w(x, y)/∂x and g(z) · ∂w(x, y)/∂y.

•
The region ahead of the notch tip where the out-of-plane shear stresses are significant progressively reduces.
In addition to this, Figure 11 makes it evident that when increasing the notch radius, the intensity of τ Max zy (z) significantly reduces as well. compared to g(z) w(x, y) / x   and g(z) w(x, y) / y   .


The region ahead of the notch tip where the out-of-plane shear stresses are significant progressively reduces. In addition to this, Figure 11 makes it evident that when increasing the notch radius, the intensity of Max zy τ (z) significantly reduces as well.  compared to g(z) w(x, y) / x   and g(z) w(x, y) / y   .


The region ahead of the notch tip where the out-of-plane shear stresses are significant progressively reduces. In addition to this, Figure 11 makes it evident that when increasing the notch radius, the intensity of Max zy τ (z) significantly reduces as well.  τ Max zy (z) as a function of the distance from the mid-plane (z/t). Different values for the notch root radius.
Subsequently, the attention was moved to discs with a thickness t = 1 mm loaded under pure mode 2 and weakened by hyperbolic notches with a root radius ρ= 0.01 mm and different notch opening angles. Results related to the plane z/t = 0.4 (where z is the distance from the mid-plane) are presented in Figures 12-14. In particular, in Figure 12, the normal stress component tangent to the notch edge, σ vv , is plotted against the polar angle θ; also, in this case, the results from the three-dimensional numerical analyses perfectly agree with the two-dimensional solution, Equations (45)-(47).  Subsequently, the attention was moved to discs with a thickness t = 1 mm loaded under pure mode 2 and weakened by hyperbolic notches with a root radius  = 0.01 mm and different notch opening angles. Results related to the plane z/t = 0.4 (where z is the distance from the mid-plane) are presented in Figures 12-14. In particular, in Figure 12, the normal stress component tangent to the notch edge, σvv, is plotted against the polar angle θ; also, in this case, the results from the three-dimensional numerical analyses perfectly agree with the two-dimensional solution, Equations (45)-(47).     In Figures 13 and 14, it is again proved that three-dimensional effects induce out-ofplane shear stresses that can be accurately predicted using the solutions derived for the antiplane deformation problem, Equations (50) and (55), independently of the material system considered. In Figures 13 and 14, it is again proved that three-dimensional effects induce out-ofplane shear stresses that can be accurately predicted using the solutions derived for the antiplane deformation problem, Equations (50) and (55), independently of the material system considered.
Finally, the disc made with material 2 and weakened by a hyperbolic notch with ρ= 0.01 mm and 2α = 90 • was loaded under pure mode III (applying the displacement fields reported in Appendix A). For this case, Figure 15 documents the existence of induced in-plane shear stresses (induced mode 2), which can be accurately predicted using Equations (45)-(47). Finally, the disc made with material 2 and weakened by a hyperbolic notch with  = 0.01 mm and 2α = 90° was loaded under pure mode III (applying the displacement fields reported in Appendix A). For this case, Figure 15 documents the existence of induced in-plane shear stresses (induced mode 2), which can be accurately predicted using Equations (45)-(47).

Conclusions
In this paper, an analytical solution for the stress fields in the close neighbourhoods of radiused notches in thick orthotropic plates under shear loading and twisting is provided. As a first step, the equations of the three-dimensional theory of elasticity were successfully reduced to two uncoupled equations in the two-dimensional space. Later, the 3D

Conclusions
In this paper, an analytical solution for the stress fields in the close neighbourhoods of radiused notches in thick orthotropic plates under shear loading and twisting is provided. As a first step, the equations of the three-dimensional theory of elasticity were successfully reduced to two uncoupled equations in the two-dimensional space. Later, the 3D stress field solution for orthotropic plates with radiused notches was presented, and its degree of accuracy was discussed by comparing the theoretical results and numerical data from 3D FE analyses.
Based on the results presented and discussed, the following main conclusions can be drawn: • Three-dimensional effects in thick plates or discs induce coupling phenomena between loading modes. In particular, out-of-plane shear stresses (mode 3) are induced on in-plane shear-loaded (mode 2) solids and can be accurately predicted using the solutions derived for the antiplane deformation problem. The vice versa is also true, independently of the orthotropic material system considered.

•
Increasing the disc or plate thickness results in an increase in the intensity of stresses induced by 3D effects; • The intensity of induced stresses is also significantly affected by the notch root radius, and the phenomenon tends to become negligible in the presence of large notch root radii.