An Approach for Stability Analysis of Gas/Fluid-Filled Cylindrical Hole in Laminated Materials

An explicit analytical workflow for cylindrical hole stability analyses in general laminated materials that possess transversely isotropic (TI) anisotropy is presented. In this approach, the calculation of the distribution of the stresses around a cylindrical hole and the failure evaluation at the hole wall consider the effects of both material elasticity anisotropy and strength anisotropy caused by material laminated structures. Material strength anisotropy is assumed to be caused by the sliding of preexisting weakness planes oriented parallel to the isotropic plane of the material. The effect of anisotropy on strength is modeled by combining a shear failure criterion for the intact matrix and a weak plane failure criterion for the planes of weakness. We derive critical pressure solutions for the stability of the intact matrix around a hole filled with gas or fluid based on the Mohr–Coulomb failure criterion and Drucker–Prager failure criterion; either one of them can be combined with the weak plane failure criterion to give the solution for hole wall shear failure pressure. The solution for hole wall fracture initiation pressure is derived based on the tensile failure criterion. This approach can be applied to holes of arbitrary orientation in general laminated materials.


Introduction
Traditional mechanical parts are mostly made of isotropic materials. Nowadays, the innovation and advancement in material production and processing methods promote the development of more sophisticated mechanical parts that can exhibit anisotropy, such as laminated composite parts, functionally graded material parts, and additive manufacturing parts. Some originally isotropic metal parts can also become anisotropic after a long period of usage due to metal fatigue or degradation. e abovementioned various types of parts often exhibit laminated features. eir elasticity can be effectively modeled as transverse isotropy, meaning that the mechanical properties are symmetric about an axis that is normal to the laminated structures. Either intrinsic or usage-induced anisotropy can make the mechanical response of a material sensitive to loading direction; thus, it is important to understand the mechanical stability of anisotropic parts before they can be used in industry applications. e mechanical concepts, theories, calculation approaches, or experimental methods established based on uniformly isotropic materials may not be directly applied for analyzing laminated materials. e inherent nonuniformity of their mechanical performance can make the mechanical analyses complicated. ere are a number of different approaches developed for studying this type of problems. Silling and Askari [1,2] introduced the dynamics theory for modeling damage and fracture in isotropic materials [1][2][3][4] and composite materials.
is approach can avoid the mathematical difficulty associated with material discontinuities by adopting an integral form of formulation in stress analysis. Chow et al. [5] presented a method of nonlinear damage analysis for anisotropic materials based on the concepts of damage surface and damage potential, and the new method was validated by its application to thin composite laminates. Sharma [6] studied the stress concentration around cutouts in an infinite layered composite plate subjected to arbitrary biaxial loading at infinity by using Muskhelisvili's complex variable method to obtain the general stress functions, from which the stress distribution can be determined. Although the complex variable method is useful in solving various elastic problems, it must be reworked when it is applied to problems with different model configurations. Zaoui et al. [7] introduced the elastic basis parameters into a Hamiltonian formulation about motion and then derived the closed-form fundamental solution of the Pastenak mathematical model based on eigenvalue analysis. Mantari et al. [8] established a new highorder shear deformation theory for the static analysis of laminated composites and sandwich plates. e cylindrical hole structure is one of the most commonly used functional structures and process structures for various types of mechanical parts. However, there is no convenient approach for determining the mechanical safety of a laminated material with holes.
Traditional numerical modeling methods, such as finite element method (FEM) and boundary element method (BEM), are commonly used in mechanical analysis of laminates. But the computational cost of three-dimensional numerical simulations can be very high when a model contains finely laminated structures. Various gradient finite element methods [9][10][11][12] have been proposed to resolve the computational efficiency issue, but these methods have their limitations and can only be applied to solve certain problems. For example, Santare et al. [9,12] benchmarked the result obtained from the gradient finite element method against the exact solution for a one-dimensional wave propagation problem and pointed out that the numerical error of finite element modeling is significant. e mechanical and mathematical modeling of laminates is still an active area of research so far.
Horgan and Chan [13] studied the stress response of a pressurized hollow cylinder made of a functionally linearly elastic isotropic material with the Young's modulus varying radially. Ying and Wang [14] proposed an analytical stress solution for a rotating multi-iron composite hollow cylinder made of radially polarized piezoelectric materials. Durodola and Attia [15] studied the deformation and stress of a functionally graded turntable and compared the modeling results obtained from a finite element method and a direct numerical integration method. Zhu et al. [16] numerically investigated the working condition of a composite sandwich cylindrical shell under hydrostatic pressure and predicted the critical failure pressure based on finite element modeling. In summary, numerical modeling of the stress distribution around a hole in a complex material is time consuming and the postanalysis for mechanical stability study is also not trivial. Considering the high complexity of model building and large computational cost of numerical modeling, practitioners would desire a fast and easy-to-use tool for stress analysis and stability evaluation for industry applications.
In this paper, we propose an analytical procedure for calculating the stresses around a cylindrical hole in general laminates and provide explicit solutions for evaluating the hole stability based on the Mohr-Coulomb failure criterion, Drucker-Prager failure criterion, weak plane failure criterion, and tensile failure criterion. e influence of the laminate internal discontinuities on mechanical stability is accounted for by incorporating the weak plane failure criterion into failure evaluation. e first two criteria (i.e., Mohr-Coulomb and Drucker-Prager) are used to model the failure of the intact matrix, and either one of them can be combined with the weak plane failure criterion to give the solution for hole wall shear failure pressure. e tensile failure criterion is used to find the fracture initiation pressure at the hole wall. ese failure criteria together can determine the safe region of the hole working pressure for given loading conditions on a part.

Methodology
ree reference coordinate systems are used in the analysis. All input parameters, such as material stiffness matrix and boundary stress tensor, are assumed to be defined in a global coordinate system (GCS), XYZ, which should align with the geometry of the material. e second coordinate system is the cylindrical hole coordinate system (HCS), in which the z-axis points downward along the hole axis and the y-axis is parallel to the isotropic plane of the material. e last coordinate system is the hole cylindrical coordinate system (CCS) that is transformed from HCS.
In this section, we will first give a brief review of the method for calculating the distribution of hole stresses in general TI media and then derive the analytical solutions for critical hole gas/fluid pressure for several commonly used failure models.

Analytical Solution for Stress Distribution around a Circular Hole in an Arbitrary TI Material.
Fang [17] presented the explicit analytical solution for stress distribution around a circular hole in an arbitrary transversely isotropic medium subjected to boundary stresses and hole internal fluid pressure. In his derivation, a hole is assumed to be infinite and homogeneous in the axial direction so that the problem can be simplified based on the generalized plain strain assumption. Hole wall is assumed to be impervious. e material surrounding the hole is assumed to possess transversely isotropy, but the elastic symmetry plane is not necessary to be parallel/normal to the hole axis, so the stiffness matrix of the material may have up to 21 nonzero components in GCS, although it only contains five independent components. Moreover, the boundary stresses may have six nonzero components because the hole and boundary principal stress directions in practices may not coincide. e constitutive relation can be expressed as where σ � [σ xx , σ yy , σ zz , σ yz , σ xz , σ xy ] T and ε � [ε xx , ε yy , ε zz , 2ε yz , 2ε xz , 2ε xy ] T are, respectively, the stress and strain vectors and the stiffness matrix C (g) that is defined in GCS has the following general form: 2 Mathematical Problems in Engineering A cylindrical hole with radius R is subjected to the hole internal gas/fluid pressure P w and the boundary stress condition that is denoted as One needs to first transform C (g) and σ (g) from GCS to HCS and then solve for the hole stress solution in HCS. Appendix gives an overview of the derivation process. In HCS, the stress solution can be written as with where the first term in equation (4), σ (0) , is the boundary stress tensor in HCS, the second term, σ (i) , which is a function of σ (g) and C (g) , is the stress alteration induced by boundary stresses, the matrix q that depends only on material property characterizes the effect of the hole gas/fluid pressure on individual stress components, E is the coordinate transformation matrix that is defined by equation (A.7). Expressions for the components in σ (i) and q are given by equations (A.40)-(A.76).
Fang [17] demonstrated that the explicit analytical expression for equation (4) can be obtained and used for hole stress calculation in arbitrary transversely isotropic media regardless of the hole orientation relative to the material elasticity symmetry axis. Moreover, this solution is valid for degenerate isotropic cases.

Solutions for Critical Hole Fluid Pressure.
In the previous section, the explicit analytical solution for stress distribution around a cylindrical hole in a general TI medium for given boundary stress conditions and hole internal gas/fluid pressure is presented. With adopted failure criteria, this hole stress solution can be used for failure evaluation of the hole wall. Explicit solutions of critical hole pressure for both shear and tensile failure are derived in this section. For shear failure, the Mohr-Coulomb and Drucker-Prager criteria are discussed separately and the weak plane approach is used to model the effect of strength anisotropy that is induced by material internal discontinuities. e stress tensor is transformed from HCS to CCS in order to perform the failure analysis: with Since the initiation of elastic failure is expected to occur at the hole wall, the calculation of critical hole gas/fluid pressure concerns the stresses at the hole wall r � R only. e components of the stress tensor σ (equation (6)) at r � R are expressed as e corresponding three principal stresses can be expressed as Note that the subscripts 1,2,3 in equations (11)-(13) do not represent their relative magnitude. σ 1 and σ 2 vary with the hole wall azimuth θ, while σ 3 is independent of θ. By inserting the expressions for σ ij into equations (11)-(13), one can obtain Mathematical Problems in Engineering e coefficients e, f, g, h, and s are expressed as

Mohr-Coulomb Shear
Failure. e Mohr-Coulomb failure criterion can be expressed as where σ min and σ max are, respectively, the minimum and maximum principle stresses, C 1 and C 2 are the material coefficients that can be related to the cohesion, S 0 , and friction angle, φ, of a material: Depending on the relative magnitude of σ 1 , σ 2 , and σ 3 , there are three possible cases: Case 1. σ min � σ 2 and σ max � σ 1 Insert equations (14) and (15) into (18) and make the expression equal to zero, which represents the failure surface; the failure criterion can be rewritten as e stable condition requires the term on the left side of equation (20) to be smaller than the other term on the right side. By taking the square of both sides of equation (20), one can obtain the following quadratic equation: with (22) e roots to equation (21) are not necessary to be the solutions to equation (20), because the square operation may introduce a fictitious root. Equation (21) may have one or two roots depending on the value of the quadratic coefficient E′. Figure 1 schematically illustrates the relationship between the expressions on the left (solid curve) and right (dashed line) sides of equation (20). e slope of the linear function on the right side of equation (20) is always negative because e < 0 and C 2 > 1. e gray region, in which the solid curve is smaller than the dashed line, represents the range of stable P w . When E′ > 0, the stable hole pressure P w is bounded within a finite range, whose upper and lower bounds are the two roots of equation (21). When E′ � 0, equation (21) reduces to a linear equation whose root defines the upper bound of safe P w . When E′ < 0, equation (20) has only one root, which is equal to the smaller root of equation (21). e larger root of equation (21) is a fictitious root that is introduced by the square operation. It corresponds to the intersection of the conjugate of the linear function on the right side of equation (20) and the other function on the left side. e lower and upper bounds of the safe P w can be written in the following concise form: P l ′ and P u ′ are functions of hole azimuth θ. e maximum P l ′ over all azimuths is the collapse pressure of hole wall, below which shear failure will occur. e minimum P u ′ over all azimuths is the allowed maximum hole pressure, above which passive shear failure will occur. e interval between the upper and lower bounds is the safe pressure region for given boundary stress conditions and material properties.
Case 2. σ min � σ 3 and σ max � σ 1 Insert equations (14) and (16) into equation (18), the failure criterion can be rewritten as Similar to the previous discussion, there are three possible situations depending on the value of the quadratic coefficient of the quadratic equation obtained by taking the square of both sides of equation (25): with (27) Figure 2 shows the relationship between the functions on the left (solid curve) and right (dashed line) sides of equation (25). Contrary to the previous case, the slope of the linear function on the right side of equation (25) is always positive because C 2 > 0 and e < 0. When E″ > 0, the two roots of equation (26) are also the roots for equation (25). When E″ ≤ 0, as shown in Figure 2, the solid curve and dashed line have only one intersection, and thus, equation (25) has one root. e lower and upper bounds of the safe hole pressure can be written in the following form: Case 3. σ min � σ 2 and σ max � σ 3 Inserting equations (15) and (16) into equation (18), the failure criterion can be rewritten as Because the slope of the linear function on the right side of equation (31) is negative, the derivation of the roots for equation (31) is similar to that for Case 1. e quadratic equation obtained through square operation is given as with Similar to the previous procedure, the lower and upper bounds of the safe hole wall pressure can be obtained as e stable region of hole pressure that fulfills the failure criterion of equation (18) can be obtained by consolidating the solutions for the above three cases. At a given hole azimuth θ, the lower and upper bounds of the stable hole pressure range, in which the hole wall will not fail along that particular azimuth, are, respectively, given as By definition, the hole critical pressure is the failure initiation pressure. e hole collapse pressure is the maximum of P (MC) l over all azimuths. Similarly, the passive shear failure pressure is the minimum of P (MC) u over all azimuths. us, the hole collapse pressure and passive shear failure pressure for the Mohr-Coulomb failure criterion can be, respectively, expressed as Equations (38) and (39) can be evaluated by simply using a grid search approach to search over the entire hole azimuth. e computational cost for finding the minimum/ maximum values is negligible, as P (MC) l (θ) and P (MC) u (θ) are given in explicit forms.

Drucker-Prager Shear Failure. Drucker-Prager shear failure criterion is written as
with where I 1 and J 2 are the stress invariants and A and B are the parameters associated with cohesion, S 0 , and friction angle, φ. Insert equations (14)- (16) into equation (40) and make the expression equal to zero; after some algebraic manipulations, equation (40) becomes ���������������������������������������� Equation (40) must have two roots so that the range of stable hole pressure is bounded. Otherwise, the solution is physically unreal. Taking the square of both sides of equation (40) and rearranging the terms, the following quadratic equation can be obtained: with When F 2 > 4EG, equation (43) has two roots that correspond to the lower and upper bounds of the safe hole pressure: (45) e hole collapse pressure and passive shear failure pressure for the Drucker-Prager failure criterion can be, respectively, expressed as Mathematical Problems in Engineering 2.5. Weak Plane Model. e previous failure models are only applicable to intact matrix. In order to model the anisotropic strength behavior of laminated materials, the weak plane approach is used to model the shear failure along the planes of weakness. e failure criterion for a plane of weakness is where τ and σ n are, respectively, the shear and normal stresses acting on the plane of weakness, S w is the intrinsic shear strength in the plane of weakness, and μ w is the sliding friction coefficient of the weak plane. It is assumed that the plane of weakness coincides with the isotropic plane of a medium. Since the y-axis is parallel to the isotropic plane (i.e., the weak plane) in HCS, the normal to the weak plane lies in the x-z plane, as shown in Figure 3. If the stiffness matrix is rotated about the y-axis in HCS by an angle θ w , which is the angle between the weak plane normal and the x-axis, the x-axis becomes normal to the isotropic plane after the rotation, resulting in the elimination of the elastic constants C 15 , C 25 , C 35 , and C 46 . Comparing the stiffness matrix before and after the rotation, one can obtain the following relationships: en, θ w can be directly determined from the elastic constants: ere are two possibilities that would lead to C 25 � C 46 � 0: (1) θ w � 0°or 90°; this indicates that the TI symmetric axis coincides with one of the coordinate axes, resulting in C 15 � C 25 � C 35 � C 46 � 0. (2) C 23 � C 12 and C 44 � C 66 ; this implies that the y-axis is normal to the isotropic plane of the formation when the formation is anisotropic, which is not possible as it contradicts the definition of the HCS. If θ w is a known parameter, then there is no need to calculate it using equation (50).
To calculate the stresses acting on the weak plane, rotate the stress tensor at the hole wall about the y-axis in HCS by an angle θ w : After the rotation, σ xx is normal to the weak plane, while σ yy and σ zz are parallel to it. e stress acting on the weak plane is given as where the normal vector n → � [1, 0, 0] T . e normal stress is e magnitude of shear stress is Inserting equations (56) and (57) into equation (48), one can obtain ����������������������������������������� P 2 w q 2 xy + q 2 xz + 2P w σ xy q xy + σ xz q xz + σ 2 Take the square of both sides of equation (58); one can obtain the following quadratic equation: with Mathematical Problems in Engineering 7 (60) Similar to the Mohr-Coulomb model, the solution for equation (58) depends on the sign of the slope of the linear term on the right side and the value of E in equation (59). Following the previous analysis procedure, one can obtain the lower and upper bounds for the safe hole pressure as for q xx ≥ 0, and for q xx < 0. e hole collapse pressure and passive shear failure pressure for the weak plane failure criterion can be, respectively, expressed as If the orientation of weak planes is different from that of the isotropic plane of the medium, one only needs to adjust the rotation matrix of equation (54) accordingly so that the stress tensor is projected onto the plane of weakness appropriately.

Tensile Fracturing Failure.
Mathematically, the tensile failure criterion can be written as where T 0 is the tensile strength that is taken to be positive. For general anisotropic materials, the magnitude of tensile strength should depend on the orientation of principal stresses relative to the anisotropy symmetric plane, as the tensile strength normal to the planes is generally different from the one parallel to them. However, the difference between tensile strength at different directions is usually small and difficult to measure, so the isotropic tensile failure criterion is used in the following derivation.
Substitute σ min with σ 2 (equation (15)) in equation (65), the tensile failure criterion becomes ������������ Note that the solution to the above equation is also constrained by the magnitude of radial stress at the hole wall: Following the same analysis procedure, the lower and upper bounds for the safe hole pressure can be obtained as z y x Hole n θ w Figure 3: Schematic illustrates the hole coordinate system (HCS), in which the z-axis is along the hole axis direction and the y-axis is parallel to the material isotropic plane that is represented by the gray lines. n → is the normal to the material isotropic plane (i.e., the material symmetry axis direction). θ w is the angle between n → and the x-axis. 8 Mathematical Problems in Engineering e minimum P (TF) u throughout the hole wall is the tensile fracture initiation pressure that is given as

Summary of the Steps for Hole Stability Analysis
(1) Divide the hole azimuth from 0 o to 180 o evenly into N intervals. e sampling points in azimuth are denoted as θ n � n · Δθ, n � 0, 1, . . . , N − 1, where Δθ is the scanning interval. Because the cylindrical hole stress variations are symmetric with respect to the hole axis, the failure evaluation only needs to be performed at azimuths from 0°to 180°. (2) Based on the given boundary stress conditions and material properties, calculate the coefficients e, f, g, h, and s (equation (17)) at θ n (n � 0, 1, . . ., N − 1). ese coefficients are independent of the hole pressure.
where P (M) L represents the matrix collapse failure pressure that is calculated from equation (38) or (45) for the Mohr-Coulomb model or Drucker-Prager model, respectively. (6) Similar to step (5), the final passive shear failure pressure can be given as where P (M) U represents the matrix passive shear failure pressure that is calculated from equation (39) or (47) for the two different shear failure models of the matrix. (7) e fracture initiation pressure is given by equation (70).

Conclusions
A complete workflow for cylindrical hole stability analysis in general laminated materials using an analytical approach is presented. Solutions of critical hole pressure for the Mohr-Coulomb model, Drucker-Prager model, weak plane model, and tensile failure model are given. is workflow can be applied to the stability analysis of arbitrary cylindrical hole in general TI media and is also applicable to degenerate isotropic media. In a TI medium, the effect of strength anisotropy is accounted for through constraining the hole collapse pressure by consolidating a shear failure model for the material matrix and the weak plane failure model. Also, it has been shown that it is necessary to consider both the lower and upper bounds of the stable hole pressure region in shear failure analysis, because collapse failure and passive shear failure can occur together at the hole wall.

Data Availability
Climate and demand data used to support the findings of this study are available from the first author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.