# Methodology for Determining the Compliance Matrix and Analyzing the Skew Anisotropic Plate with an Arbitrarily Shaped Hole

Authors: Xiaoli Zhang; Aizhong Lu; Kezhi Song; Guangming Wu; Lei Wang
Journal: Advances in Materials Science and Engineering
 (2022)

---

## Abstract

The previous analytical methods mainly focused on the analysis of orthotropic, transversely isotropic, or even isotropic materials, and little analyses have been conducted on nonorthogonal and extreme anisotropic material properties. In this study, the infinite skew anisotropic plate with an arbitrary-shaped hole are considered. Initially, two local Cartesian coordinate systems are established along the skew orientations in the plate, and the equations for the elastic constants of the material are established according to the transformation formulae under rotation of the coordinate axes. On this basis, the elastic compliance matrix for the skew anisotropic material can be obtained in combination with the engineering constants (Young’s modulus, Poisson’s ratio, and shear modulus). Subsequently, the complex variable method is utilized to derive the analytical solutions for stress along the hole boundary. The elliptical, hexagonal, and square holes perforated in skew anisotropic plates are used as examples, respectively, and the stress distributions along the hole boundary are analyzed for different external loadings. Finally, the theoretical method presented in this paper is verified by means of the Finite Element Analysis of ANSYS software.


---

## Body

## 1. Introduction

Anisotropic materials are widely used in fields such as aerospace, machinery manufacturing, civil engineering, and automobile industry owing to the excellent mechanical properties. In order to meet the requirements of assembly and structure functions, various holes are performed in anisotropic plates, which leads to a reduction in the stability and safety of the structure. It is necessary to analyze the stress distributions to ensure the sufficient strength of structure and to predict the potential risk. Therefore, many scholars at home and abroad have conducted extensive research on the perforated anisotropic structures.The complex variable method developed by Muskhelishvili [1] is an effective approach for studying the analytical solution of stress in elastic anisotropic mediums. Regarding the orthotropic plate with an elliptical hole, Lekhnitskiit [2, 3] firstly obtained accurate stress and displacement solutions for different external loadings at infinity, using the complex stress functions and conformal transformation method [4, 5]. Since then, the stress distributions in anisotropic mediums with circular [6, 7], triangular [8–10], rectangular [11–13], hexagonal [14, 15], or irregular-shaped holes [16] were investigated by scholars all over the world, subject to various external loadings. Sharma [17] even analyzed the stacking sequence of orthotropic lamina in an infinite laminated composite plate and proposed a general method to determine the stress field around a polygonal hole. Targeted for isotropic and anisotropic materials, Setiawan and Zimmerman [18] presented a unified methodology to compute the stresses around an arbitrarily shaped hole. Moreover, considering that three conformal mapping functions should be involved in solving the anisotropic problem, Manh et al. [19] derived an analytical solution for arbitrarily shaped tunnels excavated in transversely isotropic rock mass. Lu et al. [20] established three polar coordinate systems and presented an accurate analytical method to analyze the stress distribution along an arbitrarily shaped tunnel excavated in orthotropic rock mass, using the power-series method. Based on this theoretical solution, Wang et al. [21] and Zhang et al. [22] performed the optimum design of hole shape, fiber angle, and hole orientation of orthotropic plates, respectively, in order to decrease the stress concentration along the hole boundary, using a Differential Evolution (DE) algorithm. Sun et al. [23] introduced the Boundary Element Method (BEM) into the analysis of stress concentration in orthotropic media.In addition, many scholars have also carried out research on the buckling analysis or vibration response of plates with different shapes by using other solution methods. Civalek and Avcar [24] presented the free vibration and buckling analyses of functionally graded carbon nanotube-reinforced (CNTR) laminated nonrectangular plates by utilizing a Discrete Singular Convolution (DSC) method, which was also used in the investigation of free vibration response of functionally graded cylindrical shells [25]. Ng et al. [26] were concerned with the comparison of Discrete Singular Convolution and generalized Differential Quadrature (DQ) for the vibration analysis of rectangular plates. For buckling analysis of plates with different shapes, the Differential Quadrature and Harmonic Differential Quadrature (HDQ) methods were proposed by Civalek [27]. Besides, Zhu et al. [28] utilized the Finite Element Method (FEM) to study static and free vibration modeling of CNTR composites plates with the first-order shear deformation plate theory. Mishra and Barik [29] gave the nonuniform rational B-spline augmented Finite Element Method for stability analysis of arbitrary thin plates. However, the influence of hole performed in the plate has not been concerned in these researches.Moreover, although extensive theoretical studies have been carried out to derive analytical solutions for anisotropic structures, the most complex conditions involved in the abovementioned studies have mostly been orthotropic or even transversely isotropic materials. General conditions of nonorthogonal materials have not been discussed. For orthotropic materials, there are three elastic symmetry planes. The elastic compliance matrix in the generalized Hooke’s law can easily be obtained when taking the directions of the axes along the principal material directions [30]. However, coupling exists between normal stress and shearing strain, among others, for general conditions of nonorthogonal materials. Determination of the compliance matrix, especially the elements that reflect the coupling, becomes far more complicated.At present, the inverse analysis method using numerical software is utilized to identify the elastic constants for general anisotropic materials. Using the Boundary Element Method (BEM) combined with a minimization algorithm, Comino and Gallego [31] presented an inverse approach for the plane problem and determined the six elastic constants. Hematiyan et al. [32] inversely analyzed the elastic constants according to the displacement measurements from more than one static experiment and using the multiloading BEM. Regarding a three-dimensional generally anisotropic solid with arbitrary geometry, Hematiyan et al. [33] proposed a new inverse approach for the identification of all elastic constants by means of measured strain data. However, inverse analysis method is carried out by inputting the displacement or strain data measured from experiments. These methods are very susceptible to small changes in the input data [34], which makes it difficult to guarantee the accuracy of the results. Moreover, further experiments are conducted in order to address the ill-posed nature of the inverse analysis, thereby increasing the computation amount. To overcome this difficulty, it is essential to establish a simple, efficient, and accurate method to determine the elastic constants for general anisotropic materials.In this study, a methodology is provided for determining the compliance coefficient matrix of skew anisotropic materials, where only one elastic symmetry plane exists and two sets of fibers are not orthogonal but skew at any angle. On this basis, the infinite skew anisotropic plate with an arbitrarily shaped hole is considered, and the analytical stress solutions are derived using the conformal transformation method of complex variable function. The presented method is verified by comparison with the numerical simulation of ANSYS software. The influences of external loadings on stress distribution along the hole boundary are analyzed using the elliptical, hexagonal, and square holes as examples.The methodology described in this study can be applied to obtain the material constants of anisotropic plates in practical engineering. It is based on theoretical derivation and the engineering constants, which can be determined by means of simple uniaxial tensile and pure shear tests. Compared with the previous numerical and inverse analysis methods, it has the advantages of simple calculation, high precision, and universality, and it can provide a theoretical basis for future research of the anisotropic material.


## 2. Problem Description and Assumption

The only difference in the fundamental equations of anisotropic and isotropic linear elastic mechanics is the constitutive equations, because geometrical and equilibrium conditions are independent of material physical properties. For general condition of anisotropic materials, there are 13 independent elastic constants when at least one elastic symmetry plane exists, to which thez-axis of the coordinate system is perpendicular. In this case, the generalized Hooke’s law in the Cartesian coordinate system xyz is expressed as follows:(1)εx=a11σx+a12σy+a13σz+a16τxy,εy=a12σx+a22σy+a23σz+a26τxy,εz=a13σx+a23σy+a33σz+a36τxy,γyz=a44τyz+a45τxz,γxz=a45τyz+a55τxz,γxy=a16σx+a26σy+a36σz+a66τxy,where aij (i, j = 1, 2, 3, 4, 5, 6) are the elastic constants that determine the material properties. Determination of these values is a prerequisite to analyze the perforated anisotropic structures. When the material is orthotropic and the Cartesian coordinate system is established along the principal directions of the material, the number of independent elastic constants is reduced to nine. The values of a16, a26, a36, and a45 are equal to zero. The other elastic constants in (1) can be obtained in terms of engineering constants: Young’ modulus E, Poisson’s ratio ν, and shear modulus G [2]. However, several elastic constants, such as a16 and a26, are no longer equal to zero for nonorthogonal materials, regardless of the Cartesian coordinate system establishment. This study aims to determine the values of these elastic constants and obtain the compliance matrix through theoretical studies.The anisotropic material discussed here is illustrated in Figure1. The only elastic symmetry plane is the xoy-plane, to which the z-axis of the coordinate system is perpendicular. That is, the elastic properties along the z-direction are the same for each point in the material. Two sets of skew fibers exist in the xoy-plane and along the directions of x1 and x2, respectively. Similarly, each point in the material has the same elastic properties along the x1- or x2-direction, respectively. Taking the two fiber orientations as the x1- and x2-axes, respectively, the local Cartesian coordinate systems x1oy1 and x2oy2 are established. As can be seen from Figure 1, the global Cartesian coordinate system xoy overlaps x1oy1 and x2oy2 under rotation by φ1 and φ2 degrees clockwise around the z-axis, respectively. It is an orthotropic problem for φ2 − φ1 = π/2, which means that the fibers are orthogonal. Only the general situation of φ2 − φ1 ≠ π/2 is investigated in this study. In this case, neither the x1oz-plane nor the x2oz-plane is an elastic symmetry plane, that is, the skew anisotropic material, which appears frequently in the composite plates.Figure 1
Relationship between global coordinate systemxoy and local coordinate systems x1oy1 and x2oy2 in skew anisotropic materials.Based on the assumptions that the stress componentsσz, τxz, and τyz are equal to zero at each point in the plate and the other stress components are functions of x and y, the analysis of skew anisotropic plate is simplified to a plane stress problem. Thus, the generalized Hooke’s law in equation (1) can be reduced to the following equations:(2)εx=a11σx+a12σy+a16τxy,εy=a12σx+a22σy+a26τxy,γxy=a16σx+a26σy+a66τxy.In the above relations, the strain componentsγxz and γyz are equal to zero, and εz is not involved in the calculation of the other components. The elastic constants in local coordinate systems x1oy1 and x2oy2 are represented by bij and cij, respectively, and the generalized Hooke’s laws are expressed as follows:(3)εx1=b11σx1+b12σy1+b16τx1y1,εy1=b12σx1+b22σy1+b26τx1y1,γx1y1=b16σx1+b26σy1+b66τx1y1.Replacingbij, x1, and y1 with cij, x2, and y2, respectively, the other set of linear equations regarding cij can also be obtained. The elastic constants of anisotropic bodies are related to the direction of coordinates. Although aij, bij, and cij take different values, they describe the same material. As shown in Figure 1, the three Cartesian coordinate systems are related to one another by certain angles, and there is also a transformation relationship between aij, bij, and cij.


## 3. Determination of Elastic Constants

Initially, the transformation relationship between the elastic constants is analyzed under rotation of the coordinate axes. The elastic constant in the original coordinate systemxyz is denoted as dij, and the constant in the new coordinate system x′y′z′ is denoted as d′ij. dij and d′ij satisfy the following relationship:(4)d′ij=∑m=16∑n=16qimqjndmn=qimqjndmn,where i, j = 1, 2, 3, 4, 5, 6; m and n are dummy indices indicating summation. qij can be calculated according to Table 1.Table 1
Calculation rules ofqij.
1234561l211l212l213l12l13l13l11l12l112l221l222l223l23l22l23l21l22l213l231l232l233l33l32l33l31l32l3142l31l212l32l222l33l23l33l22 + l32l23l33l21 + l31l23l31l22 + l32l2152l31l112l32l122l33l13l33l12 + l32l13l33l11 + l31l13l31l12 + l32l1162l21l112l12l222l13l23l13l22 + l12l23l13l21 + l11l23l11l22 + l12l21lij = cos (xi′, xj), representing the cosine value of the angle between xi′-axis in the new coordinate system {xk′} and xj in the original coordinate system {xk}. x1′,x2′,x3′ represent x′, y′, z′. x1, x2, x3 represent x, y, z. The z-axis coincides with the z′-axis, and the new coordinate system is obtained by rotating the original coordinate system counterclockwise around the z-axis by an angle φ. Accordingly, the calculation rules of lij can be reduced to Table 2.Table 2
Calculation rules oflij.
xyzx′cosφsinφ0y′−sinφcosφ0z′001Substituting Tables1 and 2 into equation (4), the transformation formulae between elastic constants dij and d′ij can be obtained. For the skew anisotropic plate discussed in this study, bij and cij correspond to the elastic constants in the new coordinate systems, and aij corresponds to the elastic constant in the original coordinate system. Accordingly, the following equations containing aij and bij are established:(5)a11=b11cos4φ1+2b12+b66sin2φ1cos2φ1+b22sin4φ1−2b16cos2φ1+b26sin2φ1sinφ1cosφ1,a12=b11+b22−2b12−b66sin2φ1cos2φ1+b12−b26−b16cos2φ1−sin2φ1sinφ1cosφ1,a16=−2b22sin2φ1−2b11cos2φ1+2b12+b66cos2φ1−sin2φ1sinφ1cosφ1+b16cos2φ1cos2φ1−3sin2φ1+b26sin2φ13cos2φ1−sin2φ1,a22=b11sin4φ1+2b12+b66sin2φ1cos2φ1+b22cos4φ1+2b16sin2φ1+b26cos2φ1sinφ1cosφ1,a26=−2b22cos2φ1−2b11sin2φ1−2b12+b66cos2φ1−sin2φ1sinφ1cosφ1+b16sin2φ13cos2φ1−sin2φ1+b26cos2φ1cos2φ1−3sin2φ1,a66=4b11+b22−2b12−b66sin2φ1cos2φ1+b66+4b16−b26cos2φ1−sin2φ1sinφ1cosφ1.Replacingbij and φ1 with cij and φ2, respectively, the transformation formulae regarding aij and cij can also be obtained. As the above equations are not sufficient to complete the solution of the elastic constants, subsequently, equation (3) and the generalized Hooke’s law regarding bij are analyzed in the following section in combination with the engineering constants.Generally, the engineering constants are obtained by means of tests in the principal directions of the material. When a tensile test is conducted in thex1-direction, that is, only stress σx1 is applied to the structure, the first formula in equation (3) is reduced to εx1=b11σx1; thus, b11=εx1/σx1. According to the physical significance, b11 is the reciprocal of Young’s modulus Ex1 in the x1-direction; that is, b11=1/Ex1.Consideration of the second formula in equation (3) leads to εy1=b12σx1.εx1 and εy1 are substituted into the following equation:(6)νx1y1=−εy1εx1=−b12σx1b11σx1=−b12b11,where νx1y1 is Poisson’s ratio, in which the first subscript x1 indicates the applied stress direction and the second subscript y1 indicates the positive strain in the orthogonal direction.Therefore,(7)b12=−νx1y1b11=−νx1y1Ex1.If the tensile test is performed in they1-direction, only stress σy1 is applied to the structure. We obtain from the second formula in equation (3) that εy1=b22σy1, and then we have b22=εy1/σy1. Similarly, b22 is the reciprocal of Young’s modulus Ey1 in the y1-axis; that is, b22=1/Ey1.Moreover, when a pure shear test is conducted, only shear stressτx1y1 is applied. Based on the third formula in equation (3) that γx1y1=b66τx1y1, we have(8)b66=γx1y1τx1y1=1Gx1y1,where Gx1y1 is the shear modulus.In summary,bij in the above relations can be expressed by Young’s moduli, Poisson’s ratios, and shear modulus, which can be obtained by means of simple uniaxial tensile tests and a pure shear test.In the same manner, the elastic constantsc11, c12, c22, and c66 can also be obtained by analyzing the generalized Hooke’s law regarding cij, when the stresses σx2, σy2, and τx2y2 are applied, respectively.(9)c11=1Ex2,c12=−νx2y2Ex2,c22=1Ey2,c66=1Gx2y2,where Ex2 and Ey2 are Young’s moduli in the x2-axis and y2-axis, respectively; νx2y2 is Poisson’s ratio, and Gx2y2 is the shear modulus. They can also be determined by means of simple uniaxial tensile and pure shear tests.Analyzing the two sets of transformation formulae regardingaij, bij, and cij yields(10)b66−4b12=c66−4c12,b11+b22+2b12=c11+c22+2c12.Upon combining the engineering constants and the previous analysis, the unknown constants ofbij and cij are b16, b26, c16, c26, c12, and c66.Define(11)sinφ1cos3φ1=g1,sin3φ1cosφ1=g2,sin2φ1cos2φ1=g3,sinφ2cos3φ2=h1,sin3φ2cosφ2=h2,sin2φ2cos2φ2=h3.The following relations containingb16, b26, c16, c26, c12, and c66 are proposed based on equation (5).(12)−2b16g1−2b26g2+2c16h1+2c26h2−2c12h3−c66h3=c11cos4φ2−b11cos4φ1−2b12+b66g3−b22sin4φ1+c22sin4φ2,b16−b26g1−g2+c26−c16h1−h2+2h3−1c12+c66h3=c11h3−b11+b22−2b12−b66g3−b12+c22h3,b16cos4φ1−3g3+b263g3−sin4φ1−c16cos4φ2−3h3−c263h3−sin4φ2+2c12h1−h2+c66h1−h2=2c11h1+2b22g2−2b11g1+2b12+b66g1−g2−2c22h22b16g2+2b26g1−2c16h2−2c26h1−2c12h3−c66h3=c11sin4φ2−b11sin4φ1−2b12+b66g3−b22cos4φ1+c22cos4φ2,b163g3−sin4φ1+b26cos4φ1−3g3−c163h3−sin4φ2−c26cos4φ2−3h3−2c12h1−h2−c66h1−h2=2c11h2+2b22g1−2b11g2−2b12+b66g1−g2−2c22h1,4b16−b26g1−g2−4c16−c26h1−h2+8c12h3+4h3−1c66=4c11h3+c66−4b11+b22−2b12−b66g3−b66+4c22h3.Equation (12) can be expressed in the form of(13)AX=B,where(14)X=b16b26c16c26c12c66T,A=−2g1−2g22h12h2−2h3−h3g1−g2g2−g1h2−h1h1−h22h3−1h3cos4φ1−3g33g3−sin4φ13h3−cos4φ2sin4φ2−3h32h1−h2h1−h22g22g1−2h2−2h1−2h3−h33g3−sin4φ1cos4φ1−3g3sin4φ2−3h33h3−cos4φ22h2−h1h2−h14g1−g24g2−g14h2−h14h1−h28h34h3−1,B=c11cos4φ2−b11cos4φ1−2b12+b66g3−b22sin4φ1+c22sin4φ2c11h3−b11+b22−2b12−b66g3−b12+c22h32c11h1+2b22g2−2b11g1+2b12+b66g1−g2−2c22h2c11sin4φ2−b11sin4φ1−2b12+b66g3−b22cos4φ1+c22cos4φ22c11h2+2b22g1−2b11g2−2b12+b66g1−g2−2c22h14c11h3−4b11+b22−2b12−b66g3−b66+4c22h3.Providing the values of anglesφ1 and φ2 between the global and local coordinate systems, the coefficient matrix A and vector B are obtained; then, X is easy to determine. Thus far, all the bij and cij values have already been obtained. Substituting bij into equation (5) provides the elastic constants aij in the global coordinate system in equation (2). That is, the compliance matrix for skew anisotropic material is determined. It is noteworthy to mention that the proposed theoretical method is also applicable for the plane strain problem.


## 4. Derivation of Analytical Stress Solutions

Since the material constants of the skew anisotropic plate are obtained, the stress distribution around the hole can be analyzed, as presented in this section. As illustrated in Figure2, an arbitrarily shaped hole exists in the plate of Figure 1, and uniform in-plane static loadings are applied to the edge at infinity and the application of loadings is completed instantaneously. When the size of the plate is much larger than that of the hole and the hole position is not near the edge, the research domain is considered to be infinite.Figure 2
Perforated skew anisotropic plate under in-plane loadings.The conformal transformation of the complex variable function is an effective method for solving the stress problem of anisotropic plate with an irregular-shaped hole. Three complex variables (z = x + iy, z1 = x + μ1y, and z2 = x + μ2y) are involved in the solving process; thus, three polar coordinate systems (ζ = ρeiθ, ζ1 = ρ1eiθ1, and ζ2 = ρ2eiθ2) are established. With the analytic mapping functions (z = ω(ζ), z1 = ω1(ζ1), and z2 = ω2(ζ2)), the outer hole region in physical plane (z-, z1-, and z2-plane) is mapped onto the outer unit circle region in image plane (ζ-, ζ1-, and ζ2-plane), respectively. The general expressions of the mapping functions are as follows [20]:(15)z=ωζ=Rζ+∑k=0∞Ckζ−k,(16)z1=ω1ζ1=γ1Rζ1+∑k=1nCkζ1−k+δ1R¯1ζ1+∑k=1nC¯kζ1k,(17)z2=ω2ζ2=γ2Rζ2+∑k=1nCkζ2−k+δ2R¯1ζ2+∑k=1nC¯kζ2k,where Ck and R in equation (15) are coefficients denoting the hole geometry and size, respectively. Equation (15) can represent a variety of hole shapes, provided that k is sufficiently large. Herein, n is the maximum value of k.In equations (16) and (17),(18)γ1=1−iμ12,γ2=1−iμ22,δ1=1+iμ12,δ2=1+iμ22,μ1=α1+iβ1,μ2=α2+iβ2,where αk and βk (k = 1, 2) are coefficients related to material properties, and β1 > 0 and β2 > 0.The relational expressions betweenζ, ζ1, and ζ2 are established as follows:(19)γ1Rζ1+∑k=1nCkζ1−k+δ1R¯1ζ1+∑k=1nC¯kζ1k=γ1Rζ+∑k=1nCkζ−k+δ1R¯ζ¯+∑k=1nC¯kζ¯−k,(20)γ2Rζ2+∑k=1nCkζ2−k+δ2R¯1ζ2+∑k=1nC¯kζ2k=γ2Rζ+∑k=1nCkζ−k+δ2R¯ζ¯+∑k=1nC¯kζ¯−k.Equations (19) and (20) are applicable for any point in the region ⎹ζ⎸ ≥ 1, and ζ1 = ζ2 = ζ = σ = eiθ in the unit circle, which corresponds to the hole boundary.For the plane stress problem of an anisotropic plate with only one elastic symmetry plane, the body forces are neglected when the external loadings are substantially larger than the forces, and the compatibility equation for Airy’s stress functionF = F (x, y) is expressed as [30](21)a22∂4F∂x4−2a26∂4F∂x3∂y+2a12+a66∂4F∂x2∂y2−2a16∂4F∂x∂y3+a11∂4F∂y4=0,where aij are elements of the compliance matrix analyzed in Section 3.The solutions to equation (21) are connected with the roots of the following characteristic equation:(22)a11μ4−2a16μ3+2a12+a66μ2−2a26μ+a22=0.As the elastic constantsaij have been determined previously, the roots of equation (22) can be calculated. It is noteworthy that the roots are conjugate complex roots μ1, μ1¯, μ2, and μ2¯. In this study, we only investigate the situation of μ1 ≠ μ2.The stress boundary conditions are needed in the derivation process. By definingF = 2Re[F1(z1) + F2(z2)], Φ1(z1) = dF1(z1)/dz1, and Φ2(z2) = dF2(z2)/dz2, the stress boundary conditions are expressed by Φ1(z1) and Φ2(z2) as(23)Φ1z1+Φ1z1¯+Φ2z2+Φ2z2¯=f1,(24)μ1Φ1z1+μ1¯Φ1z1¯+μ2Φ2z2+μ2¯Φ2z2¯=f2.Solving for stress components is reduced to solving the two analytic functions Φ1(z1) and Φ2(z2). The research domain is infinite, and no external loadings exist along the hole boundary. Accordingly, f1 = f2 = 0, and Φ1(z1) and Φ2(z2) are written as(25)Φ1z1=B∗z1+Φ10z1,(26)Φ2z2=B′∗+iC′∗z2+Φ20z2,where B∗, B′∗, and C′∗ are real constants. By utilization of the external loadings at infinity (i.e., σx∞, σy∞, and τxy∞) and the constants α1, α2, β1, and β2, the constants can be determined as [30](27)B∗=σx∞+α22+β22σy∞+2α2τxy∞2α2−α12+β22−β12,(28)B′∗=α12−β12−2α1α2σy∞−σx∞−2α2τxy∞2α2−α12+β22−β12.(29)C′∗=α1−α2σx∞+α2α12−β12−α1α22−β22σy∞+α12−β12−α22−β22τxy∞2β2α2−α12+β22−β12.Substitution of the expressions forz1 and z2, given by equations (16) and (17), into (25) and (26), respectively, gives the following forms of Φ10z1 and Φ20z2:(30)Φ10z1=Φ10ω1ζ1=Φ1∗0ζ1=∑k=0∞akζ1−k,(31)Φ20z2=Φ20ω2ζ2=Φ2∗0ζ2=∑k=0∞bkζ2−k.Thereafter,ak and bk (k = 0, …, ∞) are determined by substituting equations (23)–(31) into the stress boundary conditions and using the power-series method.Finally, by means of equations (16), (17), and (23)–(31), we obtain(32)Φ1′z1=B∗−∑k=1nkakζ1−k−1γ1R1−∑k=1nkCkζ1−k−1+δ1R¯−1/ζ12+∑k=1nkC¯kζ1k−1,(33)Φ2′z2=B′∗+iC′∗−∑k=1nkbkζ2−k−1γ2R1−∑k=1nkCkζ2−k−1+δ2R¯−1/ζ22+∑k=1nkC¯kζ2k−1.In the Cartesian coordinate system, the stress components are determined by means of the following equations:(34)σx=2Reμ12Φ1′z1+μ22Φ2′z2,σy=2ReΦ1′z1+Φ2′z2,τxy=−2Reμ11′z1+μ2Φ2′z2.The stress componentsσρ, σθ, and τρθ in the orthogonal curvilinear coordinate system are calculated by(35)σρ+σθ=σx+σy,σθ−σρ+2iτρθ=ζ2ρ2ω′ζω′ζ¯σy−σx+2iτxy.Along the hole boundary, whereρ = 1 and ζ1 = ζ2 = ζ = σ = eiθ, the analytic solutions for stress components can easily be obtained by employing equations (32)–(35).


## 5. Examples and Comparisons

### 5.1. Comparison with Numerical Results Obtained by ANSYS Software

The Finite Element Method (FEM) is a particularly effective method for analyzing complex structural problems. In order to verify the correctness of the proposed method in this paper, an elliptical hole is used as an example, and the tangential stresses along the hole boundary are simulated by utilizing the ANSYS software in this section. Thereafter, the numerical results are compared to those obtained by means of the analytical method, using the same material parameters.For the plate in Figure2, the angles between two sets of skew fibers and the global coordinate system xoy are taken as φ1 = 10° and φ2 = 150°, respectively. The engineering constants are Ex1=60GPa, Ex2=40GPa, Ey1=15GPa, Ey2=10GPa, νx1y1=0.28, and Gx1y1=6.0GPa. On this basis, the calculation of the stress components along the hole boundary is conducted by implementing the presented theory in the FORTRAN program.The results of the material constantsaij in equation (2) are expressed as follows:(36)a11=1.426×10−2GPa−1,a12=−1.226×10−2GPa−1,a16=−1.753×10−3GPa−1,a22=8.427×10−2GPa−1,a26=4.793×10−2GPa−1,a66=0.136GPa−1.The height of the elliptical hole is 2a in the x-axis direction, while the width is 2b in the y-axis direction. Taking the values of a/b = 0.5, 1.0, and 1.5, the hole shapes are shown in Figure 3. For a/b = 1.5, the ANSYS numerical model is established as illustrated in Figure 4. The model range is 10 times the size of the width 2a, which is equivalent to the infinite domain problem in the analytical method. The values of the external loadings are taken as σx∞ = 1.0 MPa and σy∞ = 0.5 MPa, applying to the four boundaries at infinity. The four-node quadrilateral element Plane182 is adopted, and there are 252,720 nodes and 252,000 elements in the entire research domain after meshing the model by means of the mapped method. In order to improve the computational accuracy of the numerical solutions, the grid near the hole boundary is denser, while the grid away from the hole boundary is sparser. For the purpose of stress analysis, the constraints in the x- and y-directions are applied to point A, and the constraint in the x-direction is applied to point B, in order to limit the rigid body displacements of the structure.Figure 3
Shapes of the elliptical hole for different values ofa/b.Figure 4
Numerical calculation model of the ANSYS software.The tangential stresses obtained by the two methods for different ratios ofa/b are illustrated in Figure 5. The solid and dashed lines represent the analytical and numerical solutions, respectively. It is observed from Figure 5 that the results obtained by different methods are in good agreement, which verifies the theoretical method presented in this paper.Figure 5
Tangential stresses along the boundary of different elliptic holes obtained by analytical and numerical methods.


### 5.2. Stress Analysis for Different Hole Shapes and External Loadings

Elliptical, hexagonal, and square holes in the skew anisotropic plates are used as examples in this section. The values of external loadings are taken asσx∞ = σ, τxy∞ = 0, and σy∞ = λσx∞ = λσ, in which the value of σ remains unchanged and only λ is taken for different values, namely, 0.5, 1.0, and 1.5. λ is the biaxial loading factor, and λ = 1 represents the equibiaxial loading condition. The size of the hole has no influence on the calculation of the stress field, because the research domain is infinite. In order not to lose the generality, take R1 = 1.0 + 0.0i. In consideration that the hole shapes are all symmetric about the x-axis, the imaginary part of the coefficient Ck is equal to zero.The coefficients of the mapping functions for the three types of holes are given as follows [35].For the elliptical hole,a/b = 1.5, n = 1, and C1=a−b/a+b=0.2.For the hexagonal hole,n = 29, C5 = 0.0667, C11 = 0.0101, C17 = 0.0036, C23 = 0.0018, and C29 = 0.0010.For the square hole,n = 19, C3 = −0.1667, C7 = 0.0179, C11 = −0.0057, C15 = 0.0026, and C19 = −0.0014.All otherCk are equal to zero.The hole shapes illustrated in Figure6 are obtained based on the aforementioned values of Ck and R. By using the same values of skew angles and elastic constants as those in Section 5.1, the analytical solutions of the tangential stress along the hole boundary are illustrated in Figures 7–9 for different values of the biaxial loading factor, respectively.Figure 6
Different shapes of hole in the skew anisotropic plate: (a) ellipse, (b) hexagon, and (c) square.
(a)(b)(c)Figure 7
Normalized tangential stresses along boundary of the elliptical hole for different external loadings (λ = 0.5, 1.0, and 1.5).Figure 8
Normalized tangential stresses along boundary of the hexagonal hole for different external loadings (λ = 0.5, 1.0, and 1.5).Figure 9
Normalized tangential stresses along boundary of the square hole for different external loadings (λ = 0.5, 1.0, and 1.5).With no external loading applied on the hole boundary, the radial stressσρ and shear stress τρθ are equal to zero. Therefore, only the distribution of the normalized tangential stress σθ/σ, which is the ratio of the tangential stress to the external loading in the x-direction, is provided in Figures 7–9. It should be noted that angle θ in Figures 7–9 is the polar angle in the ζ-plane, which begins from the positive x-axis and rotates anticlockwise. Sign convention is defined as positive for tension and negative for compression.As can be seen from Figures7–9 that, for the perforated skew anisotropic plate, the distributions of the tangential stress exhibit asymmetry, and mainly tensile stresses are presented for the biaxial loading factor λ = 0.5, 1.0, and 1.5. Only when λ = 1.5, a small region with compressive stresses exists. However, the magnitude of the compressive stress is much smaller than that of the tensile stress. Owing to the skew anisotropic properties of the material, the region with compressive stress deviates from the intersection area of the hole boundary and y-axis (i.e., θ = 90° and 270°) and is mainly in the vicinity of θ = 80° and 260°.Moreover, when the external loading in thex-direction σx∞ remains unchanged, the tangential stresses near the intersection of the hole and x-axis gradually increase with increasing the value of σy∞. However, the tangential stresses near the intersection of the hole and y-axis decrease. It should be mentioned that there always exist four points in which the value of the tangential stress does not change with σy∞, such as point J at θ = 52.40° and point K at θ = 110.50° in Figure 7, where the stress values are always 1.2 and 2.3, respectively; the point near point D at θ = 60.95° and the point near point E at θ = 117.75° in Figure 8, where the stress values are always 7.3 and 9.0, respectively; and the point near point G at θ = 47.05° and the point near point H at θ = 131.10° in Figure 9, where the stress values are always 10.3 and 9.0, respectively. Considering that both the hole geometry and the external loadings are symmetric about the x-axis, the stress distribution at θ = 180°–360° is the repetition of that at θ = 0°–180°. Therefore, only the points with unchanged stresses at θ = 0°–180° are given in the above discussion. Accordingly, we can conclude that the external loadings in the y-direction have no influence on the tangential stresses of the four points. Furthermore, we can infer that when the external loadings are uniaxial tension stresses acting at infinity, the tangential stresses of the four points at the hole boundary are equal to zero.By comparing Figures7–9, we can also obtain that the stress distribution along the elliptical hole boundary is relatively uniform and the stress value is smaller. However, the stress concentration is very significant near the corners of the hexagonal and square holes. Moreover, for the same external loadings, a greater curvature of the hole corner leads to a larger value of the maximum tangential stress.The values of the maximum tangential stresses (σθ/σ)max and their positions θ are illustrated in Table 3 for the hexagonal and square holes under different biaxial loading factors. As can be seen from the data in Table 3, for the two types of hole, the maximum tangential stress does not exactly occur at the point of the corner but at the point near the corner. Besides, the position gets closer to the corner with increasing the value of the biaxial loading factor. Moreover, the position of the maximum stress even changes from the vicinity of point H to the vicinity of point G for the square hole under the loadings of λ = 1.5.Table 3
Values of the maximum tangential stresses (σθ/σ)max and their positions θ.
The biaxial loading factorλThe hexagonal holeThe square holeθ(σθ/σ)maxθ(σθ/σ)max0.5118.50°/298.50°10.2132.95°/312.95°14.61.0118.70°/298.70°11.3133.15°/313.15°18.21.5118.80°/298.80°12.744.95°/224.95°22.4


## 5.1. Comparison with Numerical Results Obtained by ANSYS Software

The Finite Element Method (FEM) is a particularly effective method for analyzing complex structural problems. In order to verify the correctness of the proposed method in this paper, an elliptical hole is used as an example, and the tangential stresses along the hole boundary are simulated by utilizing the ANSYS software in this section. Thereafter, the numerical results are compared to those obtained by means of the analytical method, using the same material parameters.For the plate in Figure2, the angles between two sets of skew fibers and the global coordinate system xoy are taken as φ1 = 10° and φ2 = 150°, respectively. The engineering constants are Ex1=60GPa, Ex2=40GPa, Ey1=15GPa, Ey2=10GPa, νx1y1=0.28, and Gx1y1=6.0GPa. On this basis, the calculation of the stress components along the hole boundary is conducted by implementing the presented theory in the FORTRAN program.The results of the material constantsaij in equation (2) are expressed as follows:(36)a11=1.426×10−2GPa−1,a12=−1.226×10−2GPa−1,a16=−1.753×10−3GPa−1,a22=8.427×10−2GPa−1,a26=4.793×10−2GPa−1,a66=0.136GPa−1.The height of the elliptical hole is 2a in the x-axis direction, while the width is 2b in the y-axis direction. Taking the values of a/b = 0.5, 1.0, and 1.5, the hole shapes are shown in Figure 3. For a/b = 1.5, the ANSYS numerical model is established as illustrated in Figure 4. The model range is 10 times the size of the width 2a, which is equivalent to the infinite domain problem in the analytical method. The values of the external loadings are taken as σx∞ = 1.0 MPa and σy∞ = 0.5 MPa, applying to the four boundaries at infinity. The four-node quadrilateral element Plane182 is adopted, and there are 252,720 nodes and 252,000 elements in the entire research domain after meshing the model by means of the mapped method. In order to improve the computational accuracy of the numerical solutions, the grid near the hole boundary is denser, while the grid away from the hole boundary is sparser. For the purpose of stress analysis, the constraints in the x- and y-directions are applied to point A, and the constraint in the x-direction is applied to point B, in order to limit the rigid body displacements of the structure.Figure 3
Shapes of the elliptical hole for different values ofa/b.Figure 4
Numerical calculation model of the ANSYS software.The tangential stresses obtained by the two methods for different ratios ofa/b are illustrated in Figure 5. The solid and dashed lines represent the analytical and numerical solutions, respectively. It is observed from Figure 5 that the results obtained by different methods are in good agreement, which verifies the theoretical method presented in this paper.Figure 5
Tangential stresses along the boundary of different elliptic holes obtained by analytical and numerical methods.


## 5.2. Stress Analysis for Different Hole Shapes and External Loadings

Elliptical, hexagonal, and square holes in the skew anisotropic plates are used as examples in this section. The values of external loadings are taken asσx∞ = σ, τxy∞ = 0, and σy∞ = λσx∞ = λσ, in which the value of σ remains unchanged and only λ is taken for different values, namely, 0.5, 1.0, and 1.5. λ is the biaxial loading factor, and λ = 1 represents the equibiaxial loading condition. The size of the hole has no influence on the calculation of the stress field, because the research domain is infinite. In order not to lose the generality, take R1 = 1.0 + 0.0i. In consideration that the hole shapes are all symmetric about the x-axis, the imaginary part of the coefficient Ck is equal to zero.The coefficients of the mapping functions for the three types of holes are given as follows [35].For the elliptical hole,a/b = 1.5, n = 1, and C1=a−b/a+b=0.2.For the hexagonal hole,n = 29, C5 = 0.0667, C11 = 0.0101, C17 = 0.0036, C23 = 0.0018, and C29 = 0.0010.For the square hole,n = 19, C3 = −0.1667, C7 = 0.0179, C11 = −0.0057, C15 = 0.0026, and C19 = −0.0014.All otherCk are equal to zero.The hole shapes illustrated in Figure6 are obtained based on the aforementioned values of Ck and R. By using the same values of skew angles and elastic constants as those in Section 5.1, the analytical solutions of the tangential stress along the hole boundary are illustrated in Figures 7–9 for different values of the biaxial loading factor, respectively.Figure 6
Different shapes of hole in the skew anisotropic plate: (a) ellipse, (b) hexagon, and (c) square.
(a)(b)(c)Figure 7
Normalized tangential stresses along boundary of the elliptical hole for different external loadings (λ = 0.5, 1.0, and 1.5).Figure 8
Normalized tangential stresses along boundary of the hexagonal hole for different external loadings (λ = 0.5, 1.0, and 1.5).Figure 9
Normalized tangential stresses along boundary of the square hole for different external loadings (λ = 0.5, 1.0, and 1.5).With no external loading applied on the hole boundary, the radial stressσρ and shear stress τρθ are equal to zero. Therefore, only the distribution of the normalized tangential stress σθ/σ, which is the ratio of the tangential stress to the external loading in the x-direction, is provided in Figures 7–9. It should be noted that angle θ in Figures 7–9 is the polar angle in the ζ-plane, which begins from the positive x-axis and rotates anticlockwise. Sign convention is defined as positive for tension and negative for compression.As can be seen from Figures7–9 that, for the perforated skew anisotropic plate, the distributions of the tangential stress exhibit asymmetry, and mainly tensile stresses are presented for the biaxial loading factor λ = 0.5, 1.0, and 1.5. Only when λ = 1.5, a small region with compressive stresses exists. However, the magnitude of the compressive stress is much smaller than that of the tensile stress. Owing to the skew anisotropic properties of the material, the region with compressive stress deviates from the intersection area of the hole boundary and y-axis (i.e., θ = 90° and 270°) and is mainly in the vicinity of θ = 80° and 260°.Moreover, when the external loading in thex-direction σx∞ remains unchanged, the tangential stresses near the intersection of the hole and x-axis gradually increase with increasing the value of σy∞. However, the tangential stresses near the intersection of the hole and y-axis decrease. It should be mentioned that there always exist four points in which the value of the tangential stress does not change with σy∞, such as point J at θ = 52.40° and point K at θ = 110.50° in Figure 7, where the stress values are always 1.2 and 2.3, respectively; the point near point D at θ = 60.95° and the point near point E at θ = 117.75° in Figure 8, where the stress values are always 7.3 and 9.0, respectively; and the point near point G at θ = 47.05° and the point near point H at θ = 131.10° in Figure 9, where the stress values are always 10.3 and 9.0, respectively. Considering that both the hole geometry and the external loadings are symmetric about the x-axis, the stress distribution at θ = 180°–360° is the repetition of that at θ = 0°–180°. Therefore, only the points with unchanged stresses at θ = 0°–180° are given in the above discussion. Accordingly, we can conclude that the external loadings in the y-direction have no influence on the tangential stresses of the four points. Furthermore, we can infer that when the external loadings are uniaxial tension stresses acting at infinity, the tangential stresses of the four points at the hole boundary are equal to zero.By comparing Figures7–9, we can also obtain that the stress distribution along the elliptical hole boundary is relatively uniform and the stress value is smaller. However, the stress concentration is very significant near the corners of the hexagonal and square holes. Moreover, for the same external loadings, a greater curvature of the hole corner leads to a larger value of the maximum tangential stress.The values of the maximum tangential stresses (σθ/σ)max and their positions θ are illustrated in Table 3 for the hexagonal and square holes under different biaxial loading factors. As can be seen from the data in Table 3, for the two types of hole, the maximum tangential stress does not exactly occur at the point of the corner but at the point near the corner. Besides, the position gets closer to the corner with increasing the value of the biaxial loading factor. Moreover, the position of the maximum stress even changes from the vicinity of point H to the vicinity of point G for the square hole under the loadings of λ = 1.5.Table 3
Values of the maximum tangential stresses (σθ/σ)max and their positions θ.
The biaxial loading factorλThe hexagonal holeThe square holeθ(σθ/σ)maxθ(σθ/σ)max0.5118.50°/298.50°10.2132.95°/312.95°14.61.0118.70°/298.70°11.3133.15°/313.15°18.21.5118.80°/298.80°12.744.95°/224.95°22.4


## 6. Conclusion

In this study, the skew anisotropic material, in which only one elastic symmetry plane exists and the fibers are in general nonorthogonal condition, is introduced into the problem of perforated anisotropic structure. From the view of mechanics, the stress analysis of skew anisotropic plate with an arbitrarily shaped hole is simplified to a plane stress problem in an infinite domain. A new methodology is provided for determining the compliance coefficient matrix by establishing the transformation formulae for the elastic constants under rotation of the coordinate axes and in combination with the engineering constants. On this basis, the analytical stress solutions for the skew anisotropic plate with an arbitrarily shaped hole are derived using the conformal transformation method of complex variable function.Taking the elliptical, hexagonal, and square hole as an example, the tangential stresses along the hole boundary are analyzed for different external loadings, and the following conclusions are reached.(1) The distributions of the tangential stress in the perforated skew anisotropic plate show significant asymmetry owing to the anisotropic properties of the material. (2) Mainly tensile stresses are presented under the tensile external loadings for the biaxial loading factorλ = 0.5, 1.0, and 1.5, and the region with smaller compressive stresses deviates from the intersection area of the hole boundary and y-axis. (3) When the external loadings are uniaxial tension stresses applied at infinity, there always exist four points at the hole boundary where the tangential stresses are equal to zero. (4) For the same external loadings and material parameters, a greater curvature of the hole corner leads to a larger value of the maximum tangential stress. (5) The maximum tangential stress does not exactly occur at the point of the corner but at the point near the corner; besides, the position gets closer to the corner point with the increase of the value of the biaxial loading factor.The ANSYS software is utilized to make a numerical verification of the presented method, and good agreement of stress distributions is found. The proposed solutions can provide insight into the mechanical behavior of skew anisotropic plates in practical engineering. Furthermore, they can also be applied in preliminary designs of future perforated anisotropic structures.



---

*Source: 5949398-2022-04-16.xml*