A Model accounting for Stiffness Weakening of Curvic Couplings under Various Loading Conditions

School of Energy and Power Engineering, Beihang University, Beijing 100191, China State Key Laboratory of Laser Propulsion and Application, Beijing Power Machinery Institute, Beijing 100074, China State Key Laboratory of Structural Analysis for Industrial Equipment, Department of Engineering Mechanics, Dalian University of Technology, Dalian 116023, China International Research Center for Computational Mechanics, Dalian University of Technology, Dalian 116023, China


Introduction
Rotors of modern aeroturbine engines are generally manufactured by assembling components with various materials [1]. Connection structures are introduced to connect all parts of a rotor structure system, and rotors containing these connection structures are usually defined as discontinuous rotor systems [2]. Due to the discontinuity of materials at the interface, the mechanical behavior of the connection structure is different from that of the continuous structure. In particular, an abrupt change of stiffness may exist and cause a considerable influence on the dynamic characteristics of the rotors [3][4][5][6]. erefore, in the design of a discontinuous rotor system, the stiffness weakening of the connection structures is required to be appropriately accounted for so that both the strength and dynamic requirements of the rotors under various complex working conditions are fulfilled.
Connection structures of turbine engine rotors can be classified as spline joints, bolted flange joints, and curvic couplings [7]. Extensive studies of spline joints and bolted flange joints have been carried out in the literature [8,9]. Most of them were dedicated to estimate the stiffness of joint regions [10]. Kim et al. [11] introduced four kinds of finite element models to analyse the connection stiffness of a structure with a bolted joint and made an experiment to select the best matched model. Yuan et al. [12] considered the contact interfaces of a rotor as elements with equivalent stiffness, which was determined by the Hertz contact [13] and GW model [14], and established a finite element model of a tie rod rotor. Liu et al. [15] defined two specific connection structures, i.e., a spline joint with a rabbet and a bolted flange joint with two mating cylindrical surfaces, and studied their contact states and bending stiffness employing the FEM. e stiffness with respect to geometric parameters was given, and the critical point of stiffness was found. Experimental studies were also undertaken to investigate the influences of stiffness on rotors. Zhang et al. [16] studied the effects of the bending stiffness of a splined shaft on the performance of a high-speed axial piston pump. Karpta et al. [17] proposed an experimental technique to investigate the effects of drive side pressure angle on the stiffness and measured the single-tooth stiffness of an involute spur gear experimentally. Koten et al. [18,19] conducted an experiment to verify a numerical model in the study of a combustion engine. However, due to the significant difference between curvic couplings and another two joints, the existing methods which help to determine the stiffness of the latter two connection structures cannot be fully applied to curvic couplings.
Curvic couplings, which are different from arc tooth face gears [20,21], are normally composed of convex teeth and concave teeth. Gleason Corporation has been developing the theory of the structural design of curvic couplings since 1945 [22]. With the application of Gleason Corporation grinding machines for curvic couplings, this method has been widely used in the turbomachinery industry. In this method, determining the stiffness of curvic couplings is one of the key problems for engineering design. Many researchers had been investigating methods for calculating the stiffness of curvic couplings. Bannister [23] presented design charts of the equivalent bending stiffness of curvic couplings by computing the center line slope of bending moment with the FEM. Yuan et al. [24] proposed a method calculating the position of the neutral layer of a rotor with curvic couplings. With this method, the equivalent stiffness was obtained and the influences of surface roughness, surface waviness, and preload on the stiffness and dynamic characteristics of the rotor were evaluated. Yin et al. [25] introduced a beam element to model curvic couplings, and the coefficient of the stiffness matrix was obtained by using the FEM. Gao et al. [26] introduced spring elements to express stiffness weakening of a connection structure, and the stiffness model was modified by experimental data. Gao et al. [27,28] simplified a connection structure into a mass-free hinge-bending spring to express stiffness weakening in the study of a tie rod rotor. Despite the recent progress in the study of the stiffness of curvic couplings, there are limited studies on the establishment of an analytical model, which is able to describe the relationship between the stiffness and the structure parameters of curvic couplings and to reveal the influence factors on the stiffness. e main purpose of this paper is to present an analytical approach for determining the stiffness of curvic couplings for the design of a discontinuous rotor system. Definitions and relations of the mechanical parameters in a typical curvic coupling system are firstly presented. Within this framework, a new mechanical model for analyzing the stiffness of curvic couplings is established, where a three-spring system is introduced to describe the equivalent stiffness of a pair of meshing teeth. e stiffness under various loading conditions, such as shearing, compression, bending, and torsion of curvic couplings, is derived based on the deformation compatibility of each tooth and the force equilibrium of the whole structure. Subsequently, the proposed model is adopted to determine the stiffness of curvic couplings with different structures, and the results are verified against those from the equivalent finite element analyses, showing a good match. Finally, a series of parametric studies are carried out, and the obtained results not only reveal the mechanism behind the stiffness weakening for a curvic coupling system but also demonstrate the behavior of key parameters in the design.

Definition of a Curvic Coupling System
In a typical curvic coupling of a shaft, the teeth are generally located at the end faces of annular shafts, showing a huband-spoke pattern, as shown in Figure 1(a). e shape of tooth surfaces is part of a cone, and the cross section of a tooth along the radius of the shaft is trapezoidal. Face-to-face contacts are established when meshing the convex tooth and concave tooth. e basic structure parameters of a curvic coupling are shown in Figures 1(b) and 1(c) [29], while the relationships of these parameters are listed in the last column of Table 1. More specifically, curvic couplings consist of two parts: one contains all convex teeth and the other contains all concave teeth. Each part has the same tooth number, which is even and expressed by symbol Z. e outer diameter, D, and tooth width, W, define the size of an annulus, where the teeth locate on, as shown in Figure 1(c). A point with distance S from O is defined to be the center of the pitch cone circle, whose radius is R C . e tangent from O to the pitch cone circle is defined to be the tangency radius R T , where the tooth thickness equals the space width between adjacent teeth. e number of the tooth thickness and space width enveloped by the pitch cone circle is denoted by N, which is odd. e pressure angle α, tooth width h W , and clearance c are shown in Figure 1(b). f h and f c are the user-defined coefficients for determining h W and c, respectively. e influences of the root fillet and tip chamfer of teeth are neglected here as they vary within a narrow range and have limited influence on these basic structure parameters.

A Three-Spring System-Based Model for Curvic Couplings
In order to establish the analytical approach accounting for the face-to-face contacts of curved surfaces of curvic couplings, a novel mechanical model, as shown in Figure 2(a), is employed. e model consists of an upper segment and a lower segment, both of which are considered rigid bodies. e lower segment is fixed, while the upper segment has six degrees of freedom (DOFs) in the coordinate system shown in Figure 2(a). e following methods are based on the condition that all the teeth are fully engaged under loads. e displacement of the upper segment can be expressed as 2 Mathematical Problems in Engineering where the components are the coordinates of the DOF. e contact of each meshing teeth can be simplified as a three-spring system, which is formed by two types of points and two types of springs, as shown in Figure 2(b).

Contact Points.
Contact points are defined by the midpoints of the intersections of the pitch plane and tooth surfaces, located at the circle of radius R T (see the hollow points shown in Figure 2(b)). Assuming the compression force on the Z-axis is large enough, the meshing teeth are always in contact, and thus, a pair of meshing teeth shares e coordinate of the i-th contact point is therefore expressed as

Fixed Points.
Two points are fixed on nearby surfaces of each tooth of both segments (see the solid point shown in Figure 2(b)). Each fixed point corresponds to a contact point nearby. e coordinates of the fixed points are denoted by x f− i and x f+ i , respectively, where the superscripts "− " and "+" denote the lower and upper segments, respectively.

Type I Springs.
Two springs connect the fixed points and their corresponding contact points of each tooth, as shown in Figure 3(a), describing the stiffness between the meshing teeth. When there is a displacement of the segments, the force of the type I spring on the i-th contact point can be written as where Δx f− i , Δx f+ i , and Δx c i are the displacements of the fixed points and their corresponding contact points, respectively, E − 1 and E + 1 are the stiffness of type I springs on the lower and upper segments, respectively, and d i is the normal direction of the i-th contact interface which can be written as Since small displacement of the upper segment occurs, it has negligible effect on d i . In addition, the lower segment is fixed, and thus, Δx

Type II Springs.
A curved spring connects the two contact points of a tooth, as shown in Figure 3(b), describing the coupled stiffness of two surfaces of a single tooth. When i is odd, the force of the type II spring on the i-th contact point is defined as   Figure 3: Schematics of the two types of springs: (a) type I spring; (b) type II spring.
where E − 2 and E + 2 are the stiffness of type II springs on the lower and upper segments, respectively. When i is even, the force can be written as Clearly, the stiffness of the springs is related to the structure parameters of curvic couplings and the stiffness of the corresponding tooth. e approach of obtaining the stiffness of a tooth is introduced in Section 4, where the specific equations of E − 1 , E + 1 , E − 2 , and E + 2 will be given. e following deductions are based on the assumption that these variables of the stiffness are available.
To solve equations (4), (6), and (7), the relationship between Δx f+ i and Δx c i needs to be established. Let the distances between the fixed points and their corresponding contact points tend to zeros before deformations occur, then one can obtain With the consideration that the upper segment is rigid, when an external load F is applied on the center of the upper segment, the relationship between Δx f+ i and the displacement of the upper segment, Δx 0 , satisfies where [30] and the external load F is defined as where the components F x , F y , and F z are forces and M x , M y , and M z are moments on the upper segment. ere are two groups of force balance equations [31] based on which the mechanical model is derived. One is the force balance equation for all contact points which can be given as e other is the force balance equation for the upper segment expressed as 2Z i�1 x Substituting equations (4), (6), (7), and (9) into equations (12)∼ (14), the relationship between Δx 0 and F can be deduced as where K is the stiffness matrix of curve couplings. e mathematical derivation process of the stiffness matrix K is complex, which can be found in Appendix. e components of the stiffness K can be finally expressed as follows.
Shearing component: Axial component: Bending component: Torsion component: Here, the parameters are defined as Mathematical Problems in Engineering 5 e shearing, compression, and bending stiffnesses are affected by both the type I and the type II spring stiffness, as the forces generated by the corresponding external loads are applied to the tooth surfaces on both sides. e torsion stiffness is only affected by the type I spring stiffness because there is no relative deformation between the tooth surfaces on both sides under the forces generated by the torque and the axial load due to the compatible deformation of the two surfaces.

Determination of Stiffness of a Single Tooth
e spring stiffness of a single tooth required in the model described above can be defined by the relationship between the displacements of the contact points and the distributed force on the tooth surface. Adopting the force balance on both sides of a tooth gives where E 1 and E 2 are the unknown stiffness of the two types of springs, respectively, u L and u R are the displacements at the two contact points of a tooth, respectively, and F L and F R represent the resultant force on both sides of a tooth.
As the tooth stiffness serves as a bridge connecting the structure parameters and the stiffness of curvic couplings, the relational expression between the parameters and the tooth stiffness should be obtained. In this section, the relational expression is achieved based on the stiffness analysis of a single tooth. Some basic assumptions are summarized as follows.
Assumption 1. Each tooth is discretized to many trapezoid pieces along the direction of the tooth width, as shown in Figure 4. erefore, a three-dimensional problem is converted to a two-dimensional problem.
e shearing effect between two adjacent trapezoid pieces is ignored. e deformation of each trapezoid piece in the tooth width direction is ignored. In other words, the mechanical analysis of a trapezoid piece is considered to be plane strain analysis.

Assumption 3.
Each tooth surface is subjected to distributed force, which stays unchanged in tooth height but varies in tooth width. It indicates that the force applied on each side of a trapezoid piece is uniformly distributed.

Assumption 4.
e displacement of the contact points (the point at h T still named the contact point after discretization) on the same side of trapezoid pieces of a tooth is the same.

Analysis of a Discretized Tooth
Plane. Similar to equation (21), the stiffness of a single trapezoid piece is also defined by the relationship between the displacements of the contact points and the force on the surface. According to Assumption 3, without loss of generality, the uniformly distributed forces with the sum equaling 1 are applied only on the right side with clearance c, which means that F L � 0 and F R � 1, as shown in Figure 5. en, the spring stiffness can be obtained by solving the following equation: where E W 1 and E W 2 are the stiffness of the two types of springs in one single trapezoid piece. e purpose of this section is to obtain the expressions of E W 1 and E W 2 with respect to the geometry parameters of a trapezoid piece. ere are three DOFs in the geometry of an isosceles trapezoid. However, with the same external load, F R � 1, the displacements of u L and u R are only determined by the geometry shape of the trapezoid piece, not by the specific size. erefore, only two geometry parameters are adopted in this work including the pressure angle, α, and the ratio of thickness to height, w/h W , where w denotes the thickness at the height of h T , as shown in Figure 5. e expression of w/h W is as follows: where f h is the coefficient of the whole depth h W . In engineering design, the tooth width W is usually chosen to be approximately 10% of the outer diameter D, and f h is in the range of 0.6∼1.0, as an excessively large ratio of W/D can cause grinding difficulties. e clearance c is also a parameter affecting the location of F R � 1, and its coefficient f c , which is a user-defined value, is set to be 0.12 here for simplification based on the requirement in basic engineering design.
It is difficult to determine equation (22) directly, and therefore, the data-driven approach [32] is introduced: (1) A grid sampling approach is adopted to generate n � 100 × 100 samples with different geometry parameters in the two-dimensional space, where the range of α is chosen as 20 ∘ ∼ 40 ∘ and the range of w/h W is chosen as 1.2∼2.3 according to equation (23). e samples are denoted by (2) For each 1 ≤ i ≤ n, a finite element model is established based on the geometry parameters S W i , and the software COMSOL Toolkit in MATLAB is adopted for the plane strain analysis to obtain the displacements of u L and u R . e stainless steel is adopted as the material in this work, whose Young's modulus is 2.06 × 10 11 Pa and Poisson's ratio is 0.3.
(3) en by solving equation (22) e basic interpolation models are given as where K W 1 and K W 2 are both 3 × 3 coefficient matrices. e resulted coefficients are determined as follows:    Figure 5: Schematic of the trapezoid piece.

Mathematical Problems in Engineering
can be adopted to represent E W 1 and E W 2 , as long as it has enough accuracy.

Equivalent Stiffness of a Tooth.
According to Assumption 4 that all the trapezoid pieces have the same displacement on the contact points, the resultant distributed force on the tooth surface can be calculated as where ΔR is the distance between the piece and the circle of radius R T , as shown in Figure 7(a). All the trapezoid pieces have the same values of α and h W , and the thickness on the location ΔR can be written as where "+" is for concave teeth and "− " is for convex teeth, and definitions of the two terms on the right-hand side of equation (27) can be found in Figures 7(a) and 7(b), respectively. Note that the exact expression of Δw is complex which cannot be integrated easily, and thus, an approximate form of Δw shown below is adopted, and the relative error is less than 10 − 5 compared to the exact solution of Δw.
By integrating equation (22) along the direction of the tooth width, the stiffness of a tooth can be determined. To achieve the integral, the relationship between the tooth thickness w and its corresponding radius of the circle should also be accounted for, as the tooth thickness varies with the radius. e tooth thickness w in the radius direction is equal to the arc length intercepted by the two surfaces of a tooth from the circle of the same radius, which is equal to ΔR + R T , where ΔR is the distance to the circle of radius R T . is arc length consists of two parts: one is the arc length between the two straight lines tangential to tooth surfaces through the center of the circle and the other is the arc length between the tangency line and the nearest tooth surface, as illustrated in Figure 7. e lengths are expressed by π(ΔR + R T )/Z and Δw, respectively.  (26) and (21), the final expression of spring stiffness can be obtained as where v a is a three-dimensional vector related to the structural parameters of teeth, which can be expressed as Without loss of generality, let the lower segment consist of concave teeth and the upper segment consist of convex teeth, and the superscripts "− " and "+" are used to distinguish them. erefore, the stiffness of the two springs in concave teeth is E − 1 and E − 2 , and that in convex teeth is E + 1 and E + 2 , as mentioned in Section 3.

Stiffness Analysis of Example Curvic Couplings
In this section, the proposed mechanical model is adopted to predict the stiffness of curvic couplings following the approach shown above. 12 examples of curvic couplings with For all examples, the outer diameter D is 41 mm, the teeth width W is 4.5 mm, and the tangency radius R T is 18.25 mm. Other parameters of curvic couplings, as defined in Section 2, are given in Table 2, where three levels are selected for the enveloped half pitches N, the tooth number Z, and the pressure angle α. e curvic couplings are modeled as stainless steel, which is the same material type as that included in the previous data-driven approach.

Verification of Model Predictions against FEM
Solutions. e stiffness of each example curvic coupling is obtained based on the method proposed above, and the results are listed in Table 3.
To verify the analytical solutions, 3D finite element analyses of the adopted 12 curvic coupling samples have also been performed using the commercial software ABAQUS. e material of the samples was assumed as stainless steel, with the mechanical properties being density of 7850 kg/m 3 , Young's modulus of 2.06 × 10 11 Pa, and Poisson's ratio of 0.3. e mesh consisted of 14000∼25000 8-noded hexahedral solid elements. e element size was chosen to be 0.15∼0.20 mm which was small enough to ensure the accuracy of the solutions. e DOFs of the elements below the teeth in the lower segment were set to be fully fixed, while the DOFs of the elements above the teeth in the upper segment were constrained to the origin of the local coordinate of the upper segment, as shown in Figure 8. Frictionless contacts were assumed between the teeth of the upper and lower segments for simplicity. Four finite element models have been established based on the load types to account for shearing, compression, bending, and torsion stiffnesses, respectively. e external loads, i.e., shearing, compression, bending, and torsion, were applied to the origin of the upper segment by steps, the maximum magnitudes of which were 100 N, 3000 N, 30 N·m, and 3 N·m, respectively. In the model for compression stiffness, a preload of 2000 N was applied to eliminate the influence caused by the alignment error of the elements on tooth surfaces and the stiffness was evaluated on the basis of the axial force subaccumulated from 2000 N. In other models, an axial force of 3000 N was applied to maintain the steady state of the contact under the external loads. e stiffness was determined by monitoring the displacements of the upper segment under given loads at the end of each step. e analytical solutions are compared to the FEM solutions, as shown in Table 3. e relative errors are defined by the ratios of their differences to FEM solutions. e trends of the analytical solutions are coincident with those of the FEM solutions, as illustrated in Figure 9, and the analytical solutions are less than the FEM solutions with most of the errors within 10%. e reason for the difference is that the assumptions in Section 4, such as the neglect of shear force between two adjacent trapezoid pieces and uniform distribution force on the tooth surface, which is idealized, are inconsistent with the results of the FEM.

Assessment of Stiffness Weakening.
According to the theory of material mechanics [33], the shearing, bending, compression, and torsion stiffnesses of continuous structures can be approximated by the following equations, respectively: e stiffness weakening caused by the connection of curvic couplings is evaluated by the ratio of the stiffness difference between curvic couplings and continuous structures to the stiffness of continuous structures with the same size. e results are shown in Table 4. e stiffness weakening of shearing, compression, bending, and torsion of curvic couplings is the ratio of 50.90%∼71.54%, 53.78%∼ 93.72%, 69.74%∼93.93%, and 35.60%∼44.13% to the stiffness of continuous segments, respectively. e weakening trend of bending stiffness is consistent with that of compression stiffness.

Parametric Studies
is section studies the influence of structural parameters of curvic couplings on the connection stiffness, in order to provide an insight into the optimization of the stiffness in engineering design.

Effects of the Structure Parameters on Tooth Stiffness.
According to the analysis in Section 4, the stiffness of a trapezoid piece is only determined by the geometric shape of the trapezoid piece, not by the specific size, which infers that the tooth number of curvic couplings has no effect on the stiffness of a single tooth. After eliminating the influence of the tooth number, the results from Table 3 are used to plot stiffness curves, as shown in Figure 10, where the enveloped half pitches have little influence on the tooth stiffness. Additional examples with pressure angles, i.e., 22.5°, 25°, 27.5°, 32.5°, 35°, and 37.5°, are designed according to the design theory of curvic couplings [29]. As the influence of enveloped half pitches on the stiffness of a single tooth could be ignored, the tooth number and enveloped half pitches of the examples, Z � 24 and N � 21, are selected. e results are shown in Figure 11. e differences in stiffness between convex and concave teeth are exceedingly small, while the pressure angle has a significant influence on the tooth stiffnesses E 1 and E 2 irrespective of other variations. e bending and shearing deformations of a trapezoid piece, which are caused by the transverse component of the applied unit force, account for a high proportion in the total value of u R . It should be noted that u R is inversely proportional to E 1 . With an increase in the pressure angle, the transverse component of the applied unit force decreases, which consequently decreases u R and then increases E 1 . e variation of E 2 is related to the relative deformation of the tooth surfaces on both sides, which is less affected by the increase in the pressure angle, as the compatible deformation of the two tooth surfaces leads to less relative deformation under the applied unit force on one side of the tooth. erefore, it can be inferred that the pressure angle is the main factor affecting the stiffness of a single tooth, and E 1 is more sensitive to the pressure angle than E 2 .

Effects of the Structure Parameters on the Stiffness of Curvic
Couplings. According to the design theory of curvic couplings, the size of the trapezoid pieces is inversely proportional to the tooth number when other parameters stay unchanged. e increase of tooth number implies the increase of the spring number in the spring system, which causes the increase of stiffness of curvic couplings for these springs in parallel. e stiffness of a curvic coupling is chosen to be the basic value, and the relative variation of stiffness is obtained by calculating the ratio of other stiffnesses to it. Adopting the same pressure angle and choosing the stiffness of the curvic coupling with 12 teeth as the basic value, the stiffness of curvic couplings increases almost linearly with the increase of the tooth number, as shown in Figure 12. According to the definition of curvic couplings in Section 2, when the outer diameter, tangency radius, and tooth width stay unchanged, increasing the number of teeth results in a smaller tooth size. As stated in Section 4.1, the displacements at the two contact points of a tooth only depend on the geometry shape of the trapezoid piece, not the specific size under the action of unit force. erefore, with an increasing tooth number, the stiffness of curvic couplings also increases. Additional examples with the same pressure angle as in Figure 11 are designed with the tooth number Z � 16 and enveloped half pitches N � 13. e results are shown in Figure 13. Given the same tooth number and choosing the stiffness of the curvic coupling with pressure angle 20°as the basic value, the stiffness of curvic couplings increases with the increase of pressure angle, and the compression and bending stiffnesses are more sensitive than the shearing and torsion stiffnesses, as shown in Figure 13.  1  2  3  4  5  6  7  8  9  10  11  12  Tooth number  Z  12  12  12  16  16  16  20  20  20  24  24  24  Enveloped half pitches  N  7  9  11  11  13  15  15  17  19  19  21  23  Pressure angle ( o )  α  20  30  40  30  40  20  40  20  30  20  30     Note. e unit for shearing and compression stiffnesses is 10 3 N/mm, while that for bending and torsion stiffnesses is 10 6 N·mm/rad.

Conclusions
is paper focuses on establishing the relationships between the structure parameters and the stiffness of curvic couplings and presents an analytical approach to simulating the stiffness of curvic couplings in the design of an aeroturbine engine rotor under complex loading conditions, e.g., shearing, compression, bending, and torsion.
e key conclusions can be summarized as follows: (1) A three-spring system consisting of two types of springs is proposed to describe the equivalent stiffness of a pair of meshing teeth, which is determined by the stiffness analysis of the discretized tooth plane. (2) e stiffness of a single tooth increases with the increase of the pressure angle, while the enveloped half pitches and tooth number have little influence on it, which means the influence of the latter two on the stiffness of a single tooth could be ignored. (3) e stiffness of curvic couplings is significantly smaller compared to that of continuous structures, and the pressure angle and tooth number are the major factors affecting the stiffness of curvic couplings, while the compression and bending stiffnesses are more sensitive to the pressure angle than other stiffnesses. (4) e curvic couplings with more teeth could obtain greater connection stiffness when keeping other parameters unchanged, which provides a way of optimizing the connection stiffness without changing the diameter of the shaft of an aeroengine rotor provided the parameters changed still satisfy the grinding conditions of curvic couplings. (5) e stiffness weakening of curvic couplings is tremendous, which may have enormous influence on the dynamic characteristics of the rotors. e analytical method proposed is of guiding significance to the structural design of the curvic couplings in engineering design of aeroengine rotors.
Appendix e mathematical deduction process of the stiffness matrix K in equation (15) is introduced in the following text.
By substituting equation (4), (6), (7), and (9) into equation (12), one can obtain Define C as the displacements of all the contact points on d i : and define matrix D as  According to the definition of the matrices D and W − , the balance equations of equations (13) and (14) can be rewritten in the vector form as To express K, the matrix V is defined as (A.13) Several properties of V are given as follows. e columns of V are orthogonal vectors, and the following equation holds: (A.14) e relationship between V and D is where T 0 is defined as (A. 16) In addition, W + V and W − V can be calculated as where where k + , k − , k 0+ , and k 0− are all 2 × 2 matrices given by k + � e + − e + 2 cos(π/Z) e + 2 sin(π/Z) e + 2 sin(π/Z) e + + e + 2 cos(π/Z) According to the definition of T − and T + , it can be calculated that

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

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