Abstract

Curvic couplings are frequently used in aeroengine rotors. The stiffness of the curvic couplings is of guiding significance to the engineering design of aeroengine rotors as it is significantly different from that of continuous structures. In this paper, definitions and relations of the structure parameters for a curvic coupling are firstly introduced. Based on this proposed mechanical framework, a novel mechanical model accounting for the stiffness weakening under shearing, compression, bending, and torsion is developed for curvic couplings. In this model, a three-spring system, which consists of two types of springs, is adopted to describe the equivalent stiffness of a pair of meshing teeth of curvic couplings. The spring stiffness is obtained by employing the plane strain analysis of a discretized tooth with trapezoid pieces. Subsequently, the stiffness matrix of curvic couplings is deduced based on the deformation compatibility of each tooth and the force balance of the whole structure. A series of analyses of curvic couplings with various structure types are performed to demonstrate the mechanism behind the proposed model, and the results are verified against those obtained from finite element analyses. It is shown in this study that the pressure angle is the major factor affecting the stiffness of curvic couplings, while the compression stiffness and bending stiffness are more sensitive than other stiffnesses. Furthermore, the stiffness of curvic couplings is considerably smaller compared to that of continuous structures, indicating the importance of appropriate modelling of stiffness weakening in the design of aeroengine rotors.

1. 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 [36]. Therefore, 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. The 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.

The 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. The 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.

2. 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 hub-and-spoke pattern, as shown in Figure 1(a). The 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.

The 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. The 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 RC. The tangent from O to the pitch cone circle is defined to be the tangency radius RT, where the tooth thickness equals the space width between adjacent teeth. The number of the tooth thickness and space width enveloped by the pitch cone circle is denoted by N, which is odd. The pressure angle α, tooth width hW, and clearance c are shown in Figure 1(b). fh and fc are the user-defined coefficients for determining hW and c, respectively. The 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.

3. 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. The model consists of an upper segment and a lower segment, both of which are considered rigid bodies. The lower segment is fixed, while the upper segment has six degrees of freedom (DOFs) in the coordinate system shown in Figure 2(a). The following methods are based on the condition that all the teeth are fully engaged under loads. The displacement of the upper segment can be expressed aswhere the components are the coordinates of the DOF.

The 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).

3.1. 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 (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 the same contact point. The angular position of the i-th contact point can be defined as

The coordinate of the i-th contact point is therefore expressed as

3.2. 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. The coordinates of the fixed points are denoted by and , respectively, where the superscripts “−” and “+” denote the lower and upper segments, respectively.

3.3. 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 aswhere , , and are the displacements of the fixed points and their corresponding contact points, respectively, and are the stiffness of type I springs on the lower and upper segments, respectively, and 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 . In addition, the lower segment is fixed, and thus, .

3.4. 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 is odd, the force of the type II spring on the i-th contact point is defined aswhere and are the stiffness of type II springs on the lower and upper segments, respectively. When 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. The approach of obtaining the stiffness of a tooth is introduced in Section 4, where the specific equations of , , , and will be given. The 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 and 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 obtainWith 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 and the displacement of the upper segment, , satisfieswhere [30]and the external load F is defined aswhere the components Fx, Fy, and Fz are forces and Mx, My, and Mz are moments on the upper segment.

There 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

The other is the force balance equation for the upper segment expressed as

Substituting equations (4), (6), (7), and (9) into equations (12)∼(14), the relationship between and F can be deduced aswhere K is the stiffness matrix of curve couplings. The mathematical derivation process of the stiffness matrix K is complex, which can be found in Appendix. The 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

The 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. The 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.

4. Determination of Stiffness of a Single Tooth

The 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 giveswhere E1 and E2 are the unknown stiffness of the two types of springs, respectively, and are the displacements at the two contact points of a tooth, respectively, and and 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. Therefore, a three-dimensional problem is converted to a two-dimensional problem.

Assumption 2. The shearing effect between two adjacent trapezoid pieces is ignored. The 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. The displacement of the contact points (the point at hT still named the contact point after discretization) on the same side of trapezoid pieces of a tooth is the same.

4.1. 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  = 0 and , as shown in Figure 5. Then, the spring stiffness can be obtained by solving the following equation:where and are the stiffness of the two types of springs in one single trapezoid piece.

The purpose of this section is to obtain the expressions of and with respect to the geometry parameters of a trapezoid piece. There are three DOFs in the geometry of an isosceles trapezoid. However, with the same external load, , the displacements of and are only determined by the geometry shape of the trapezoid piece, not by the specific size. Therefore, only two geometry parameters are adopted in this work including the pressure angle, α, and the ratio of thickness to height, , where denotes the thickness at the height of , as shown in Figure 5. The expression of is as follows:where fh is the coefficient of the whole depth hW.

In engineering design, the tooth width W is usually chosen to be approximately 10% of the outer diameter D, and fh is in the range of 0.6∼1.0, as an excessively large ratio of W/D can cause grinding difficulties. The clearance c is also a parameter affecting the location of , and its coefficient fc, 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 and the range of is chosen as 1.2∼2.3 according to equation (23). The samples are denoted by .(2)For each , a finite element model is established based on the geometry parameters , and the software COMSOL Toolkit in MATLAB is adopted for the plane strain analysis to obtain the displacements of and . The stainless steel is adopted as the material in this work, whose Young’s modulus is 2.06 × 1011 Pa and Poisson’s ratio is 0.3.(3)Then by solving equation (22), the responses of the spring stiffness are calculated based on the obtained displacements, which are denoted as and , as shown in Figure 6.(4)With the training data of , , and , the mathematical expressions of and with respect to the geometry parameters are obtained by least-squares interpolation. The basic interpolation models are given aswhere and are both 3 × 3 coefficient matrices. The resulted coefficients are determined as follows:

The maximum relative errors of and are 0.69% and 3.43%, respectively, while the average relative errors are 0.18% and 0.71%, respectively.

The sensitivities of the stiffness of shearing, compression, bending, and torsion to the stiffness of thetwo types of spring are analysed. With an increase of by 5%, the stiffnesses increase by 2.44%, 1.41%, 1.42%, and 2.44%, respectively. However, the stiffnesses only increase by 0.002%, 1.05%, 1.05%, and 0%, respectively, when increases by 5%, showing weak sensitivities. An alternate interpolation model can be adopted to represent and , as long as it has enough accuracy.

4.2. 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 aswhere is the distance between the piece and the circle of radius , as shown in Figure 7(a). All the trapezoid pieces have the same values of and , and the thickness on the location can be written aswhere “+” 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 is complex which cannot be integrated easily, and thus, an approximate form of shown below is adopted, and the relative error is less than 10−5 compared to the exact solution of .

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 and its corresponding radius of the circle should also be accounted for, as the tooth thickness varies with the radius.

The tooth thickness 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 , where is the distance to the circle of radius . This 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. The lengths are expressed by and , respectively.

By comparing equations (26) and (21), the final expression of spring stiffness can be obtained aswhere 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. Therefore, the stiffness of the two springs in concave teeth is and , and that in convex teeth is and , as mentioned in Section 3.

5. 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 different structural types have been selected in this study. For all examples, the outer diameter D is 41 mm, the teeth width W is 4.5 mm, and the tangency radius 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 α. The curvic couplings are modeled as stainless steel, which is the same material type as that included in the previous data-driven approach.

5.1. Verification of Model Predictions against FEM Solutions

The 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. The material of the samples was assumed as stainless steel, with the mechanical properties being density of 7850 kg/m3, Young’s modulus of 2.06 × 1011 Pa, and Poisson’s ratio of 0.3. The mesh consisted of 14000∼25000 8-noded hexahedral solid elements. The element size was chosen to be 0.15∼0.20 mm which was small enough to ensure the accuracy of the solutions. The 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. The 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. The stiffness was determined by monitoring the displacements of the upper segment under given loads at the end of each step.

The analytical solutions are compared to the FEM solutions, as shown in Table 3. The relative errors are defined by the ratios of their differences to FEM solutions. The 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%. The 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.

5.2. 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:

The 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. The results are shown in Table 4. The 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. The weakening trend of bending stiffness is consistent with that of compression stiffness.

6. Parametric Studies

This 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.

6.1. 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. The results are shown in Figure 11. The differences in stiffness between convex and concave teeth are exceedingly small, while the pressure angle has a significant influence on the tooth stiffnesses E1 and E2 irrespective of other variations. The 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 . It should be noted that is inversely proportional to E1. With an increase in the pressure angle, the transverse component of the applied unit force decreases, which consequently decreases and then increases E1. The variation of E2 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. Therefore, it can be inferred that the pressure angle is the main factor affecting the stiffness of a single tooth, and E1 is more sensitive to the pressure angle than E2.

6.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. The 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. The 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. Therefore, 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. The 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.

7. Conclusions

This 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. The 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)The 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)The 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)The 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)The stiffness weakening of curvic couplings is tremendous, which may have enormous influence on the dynamic characteristics of the rotors. The analytical method proposed is of guiding significance to the structural design of the curvic couplings in engineering design of aeroengine rotors.

Appendix

The 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 :and define matrix D aswhere

Equation (A.1) can be rewritten in the vector form aswhere and are column transformation matrices of two adjacent columns beginning with the first and second columns, respectively. For example, the 6-order column transformation matrices are written as

Define matrix

Based on equation (A.24), one can obtain

According to the definition of the matrices D and , the balance equations of equations (13) and (14) can be rewritten in the vector form as

By substituting equation (A.8) into (A.9), one can obtainwhere

To express K, the matrix V is defined aswhere

Several properties of V are given as follows. The columns of V are orthogonal vectors, and the following equation holds:

The relationship between V and D iswhere is defined as

In addition, and can be calculated aswherewhere , , , and are all 2 × 2 matrices given by

Based on equation (A.17), it can be deduced thatand then

By substituting equations (A.15) and (A.17)∼(A.18) into equation (A.11), the stiffness matrix can be deduced as

According to the definition of and , it can be calculated thatwhere

By substituting equations (A.23) and (A.16) into equation (A.22), the expression of the stiffness matrix K can be obtained, as expressed in equations (16)∼(19).

Data Availability

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

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to thank the financial support by the Ministry of Industry and Information Technology of China (JSZL2016204B102), the National Natural Science Foundation of China (51575022, 11432003, 11602050, and 11702053), and the 111 Project of China (B14013).