Application of Isogeometric Analysis Method in Three-Dimensional Gear Contact Analysis

1  Introduction

The gear contact analysis is a complicated nonlinear problem. The classical Hertz formula, first formulated by Hertz [1] in 1881, utilized mathematical and elastic methods to address a variety of contact issues associated with basic shapes and uniform contact surfaces. However, several conditions must be met for the classical Hertz contact method to address certain engineering challenges. Firstly, the traditional finite element analysis (FEA) method for the gear contact problem is to use the C(0) continuous Lagrangian basis function for boundary fitting. This often results in discontinuities at the fitted surface boundary and a failure to capture the contact surface’s exact shape. Such limitations hinder improvements in contact stress calculation accuracy. Secondly, in the pre-processing stage of contact analysis, the traditional finite element method uses tetrahedron or hexahedron meshes for mesh generation and approximates the original model with an approximate geometric mesh model at the contact boundary. These differences often reduce analytical efficiency and can result in large variances in the results, especially when analyzing the contact problems of intricate surfaces.

In 2004, isogeometric analysis (IGA) was first proposed by Hughes et al. [2] as an analysis method based on the Non-Uniform Rational B-Splines (NURBS) parametric model to unify the design model and the analysis model and then achieve a seamless connection between design and analysis. The model for IGA is composed of NURBS surfaces or volumes, so the entire model remains smooth, avoiding the problem of inaccurate model representation produced by the standard finite elements’ discrete computational unit densities. By leveraging the NURBS basis, the IGA method avoids the error of the FEA method from the geometric source.

However, despite nearly two decades elapsing since IGA’s inception, its application to complex mechanical components like gears remains limited. Several reasons have contributed to this limitation. One reason is that the complex mechanical parts with many engineering features are very hard to construct to be suitable for IGA. Standard mechanical components, including gears, bearings, bolts, flat keys, splines, and springs, are characterized by specific curves or surfaces. Representing these components as volume parametric solids is far from straightforward. The other reason is that IGA cannot yet work on the analysis tasks under complex work conditions. As a new and promising analysis method to replace the current FEA method, IGA has its gaps. What is more, in using the existing methods of analysis methods some analysis problems remain unsolved. Just as with the problem of gear contact analysis, although there are many analysis methods, including the FEA methods, finding a precise, efficient, economical, and fast analysis method remains to be solved and has attracted many researchers’ interest for decades. The basic idea of IGA is that of using the same mathematical language to express geometric shapes and physical fields. Using a planar NURBS surface or a three-variable NURBS parameter volume as the computational domain, wherein the computational units are the surface or volume units corresponding to the node interval, the NURBS basis function is used as the shape function, and the physical properties of the control vertices are used as the invariant. Therefore, for 3D models with complex structures, the key to successful geometric analysis is to decipher how to build a volume-parameterized model that meets the continuity requirements. To achieve this, it is necessary to solve two problems: the NURBS multi-piece (block) splicing method and boundary continuity coupling. Here, we used IGA to establish the three-dimensional contact model of the gear and compared its tooth profile with that of the standard involute cylindrical gear for geometric error analysis. Then, we developed a new geometric contact analysis algorithm to deal with the frictionless contact problem of 3D gears and compared the analysis results with the results obtained from theoretical calculations and commercial software.

The paper is outlined as follows: In Section 2, an overview of gear contact analysis, isogeometric analysis, and isogeometric contact analysis is presented. In Section 3, a common isogeometric contact analysis for three-dimensional calculations is derived and realized. To study the correctness of calculation formulas and processes, a standard one-fourth semicircle contact model is established, and the results of the isogeometric contact analysis are compared with the Hertz theoretical values and the empirical values from commercial software. In Section 4, a NURBS representation of a gear contact model is analyzed, and the contact Von-Mises stress is obtained and compared with the corresponding values obtained from some commercial software. Section 5 presents the conclusions and ends this paper.

The contributions of our work can be summarized as:

1) The calculation formulas for common 3D isogeometric contact analysis were derived.

2) Based on the requirements for smoothness and continuity of gear contact boundaries, a 3D gear body parameterized model suitable for geometric analysis was constructed using the NURBS multi-block splicing and Nitsche’s method. The correctness and robustness of this method were verified by comparisons with FEA results.

3) The advantages of IGA methods in 3D contact analysis were explored, including the finding that under the same calculation accuracy, there were fewer calculation units and degrees of freedom, smoother expression of contact surfaces, smoother contact stress, and significantly reduced non-physical oscillations.

2  Related Works

2.1 Contact Analysis of Gears

Gears are some of the most used parts in transmission equipment, so the gear contact problem has gained many researchers’ interest, including other problems such as the fatigue life analysis of gears [3], transmission error analysis of gears [4], and stress analysis of gears [5]. The FEA method has been widely used for the analysis of gears. Li [6] performed a finite element frictionless contact analysis on a three-dimensional gear. Mao [7] used the micro-geometric correction method to reduce the fatigue wear of gear tooth surfaces and calculated the contact stress of gear tooth surfaces with the finite element method. Gonzalez-Perez et al. [8] proposed a method to reduce the calculation cost of the gear meshing by using multi-point constraints at different sections of the gear geometry. Tamarozzi et al. [9] analyzed two gears meshing in an ad-hoc flexible multibody model and proposed an online selection of residual attachment modes for accurate local deformation prediction. Chang et al. [10] calculated the overall deformation of the tooth body by the analytical method and calculated the local deformation in the contact area by the finite element method, which greatly shortened the calculation time. Li [11,12] studied the meshing stiffness and stresses of the straight cylindrical gears by using the finite element method, which considered tooth profile modification, manufacturing errors, and other factors. Wang [13] proposed a method that could correct the tooth profile by using tooth surface contact and established a dynamic model of helical gears. Litvin et al. [14,15] used a gear geometric contact analysis technology to modify the tooth surface of helical gears in three dimensions and reduce the amplitude of transmission error, and also developed a six-degree-of-freedom helical gear dynamic model to study the time-varying meshing stiffness (TVMS) calculation method by considering axial deformation. Lin et al. [16] proposed a transmission error analysis method for helical gear systems that considered machining assembly and correction errors, and established a bent-twist shaft coupling dynamic model for the analysis and control of gear system vibration and noise.

2.2 Application of IGA

Since Hughes et al. [17] proposed the IGA method, many researchers have further developed it. Currently, this method has been applied in many fields, such as heat conduction problems [18], fluid flow problems [19,20], acoustic problems [21], electromagnetic problems [22], and fluid structure coupling problems [23–28]. IGA has also been applied in the medical domain [29–30].

IGA has been deeply and widely verified in the field of mechanical engineering. The most obvious advantage of the IGA method is its geometric accuracy, which can accurately express the geometric model. By bonding complex models with multiple NURBS surfaces, Kiendl et al. [31,32] solved the continuity problem between multiple NURBS surfaces. Based on the advantages of IGA, it has been applied to plate and shell analysis by many researchers. To study the vibration of thin-walled beams, Carminelli et al. [33] used B-splines as elements in establishing the thin-walled shell elements of multi-layer composite materials. Uhm et al. [34] used T-splines to establish a plate and shell model for IGA, and effectively solved the problem of the degree of freedom that could not be rotated at the control points. Hirschler et al. [35] proposed a dual-domain decomposition algorithm to accurately analyze the unqualified polyhedral Kirchhoff-Love shells. Jamires et al. [36] proposed an IGA-based method to study the stability of laminates and shallow shells. The proposed method was applied to the stability analysis of various examples and obtained excellent results. Norouzzadeh et al. [37] applied the IGA method to the dynamic analysis of the shell and used the method to generate the geometry of the surface area and the approximate displacement field in the shell.

The IGA method has also been combined with other research methods to solve broader problems. Chen et al. [38] realized a parametric optimization for complex structural shapes by using IGA and verified the effectiveness of this method with several complex planar shapes. Xia et al. [39] proposed a coupling model based on IGA and bond-ed peridynamics (PD) (IGA-PD) and proved that in the IGA-PD coupling model, the PD node could be constructed by the control network of IGA, which improved the calculation accuracy and efficiency of crack problems and provided a solution method for crack problems with an accurate geometric model. Aung et al. [40] combined shape optimization with IGA and applied this method to the analysis and optimization of welds. The analysis and optimization results showed that the stress concentration of welds was greatly reduced by using this method.

Since the gear is a typical part in mechanical engineering, here, we give a detailed overview. Yusuf et al. established a three-dimensional gear model with NURBS to study the strength and deformation of a gear tooth under load. By comparing the results of IGA with the corresponding values of the Lewis equation and FEA, it was proved that the IGA model could obtain more accurate gear analysis results. However, the authors only gave the result of stress analysis of a single tooth, instead of the teeth that kept contact. In addition, the established IGA gear model was not compared with the standard gear model for the geometric error of the tooth profile. As there is no specific derivation formula for the analysis of isogeometric gears [41]. Considering the mesh interference and profile modification under loads, Beinstingel et al. used NURBS to design an accurate tooth profile to explore the influence of mesh density on gear meshing stiffness matrix, then proved that the analysis established by IGA had good applicability and stability. However, the influence of the stiffness matrix derived under the quasi-static state on the gear force was not discussed [42].

For the contact analysis between mechanical parts, some researchers have also conducted relevant research with IGA. To quickly search for the control points on isogeometric models, Stadler et al. [43] proposed an efficient search algorithm for adaptive fine mesh on adjacent elements search. Temizer et al. [44] used the NURBS contact discretization node-to-surface (KTS) algorithm to study frictionless contact problems. It qualitatively showed that the application of this algorithm could improve the analysis accuracy and convergence, which focused on the effect of the proposed new algorithm on mesh discretization without analyzing its stress results. Lu [45] was the first to develop the IGA algorithm framework for contact problems. Two contact algorithms were used to analyze Hertz and interference contact examples, and the analytical solutions were compared. Mathisen et al. [46] simulated the large deformation contact problem between pipeline and trawl gear based on the IGA of mortar, and the numerical example demonstrated that isogeometric surface discretization delivered more accurate and robust predictions of the response compared to Lagrange discretization. Based on the local maximum entropy (LME) approximation, Greco et al. [47] explored the simulation of small deformation contact problems with a coupled isogeometric meshless formulation. Kamensky et al. [48] used IGA to derive the contact equation based on volume potential energy and used the structural mechanics of the atrioventricular valve to illustrate the validity of this equation. Emad et al. [49] proposed an adaptive method to deal with contact problems. The effectiveness of this method was verified by the Hertz problem and three hyperelastic contact problems, and the results showed that IGA had good convergence not only for contact pressure, but also for the limit of contact area.

In summary, although isogeometric analysis has obvious advantages in geometric accuracy, more work is needed to effectively utilize IGA in solving the problems of gears:

1) There is no specific mechanical formula for the meshing contact of two cylindrical gears combined with the theory of isogeometry.

2) Compared with the FEA method, the advantages of the IGA method in modeling and analyzing two cylindrical gears are not clear.

3  Realization of Isogeometric Contact Analysis

The process of IGA is similar to FEA, but there are also differences. For instance, while FEA uses the Lagrangian Function as its basis function, and the boundary of the model can only reach C(0) continuous, IGA uses the NURBS basis function which can reach at least C(1) continuous. The basic process of isogeometric contact analysis is shown in Fig. 1. In the whole process, definitions of the contact pairs and constraint boundary conditions are the most important components. Accordingly, the following subsections detail the contact area search, defines the contact pairs, and presents the constraints.

Figure 1: Process of isogeometric contact analysis

3.1 Search of the Contact Areas

Since the isogeometric model was presented with NURBS, searching the contact areas meant searching the control points. So, the searching operation was divided into two steps. The first step, called a global search, was to determine which elements would be contacted, and the second step called local search was to determine the contact control points on the contacted units. The difficulty lay in the local search step. As shown in Fig. 2, the spatial distance function S required the shortest distance from the spatial point to the surface.

Figure 2: Distance from a space point to a surface

As shown in Fig. 2, given a point Q that was not on the NURBS surface, its vertical foot point P0 could be found by an iterative method. Firstly, any point P1(u,v) on the surface was taken as optional and its two partial derivatives could be obtained as ru, rv. The direction vector S was given as:

S=Arω+Bru+Crv(1)

where, rw is the outer normal direction of point P1 on the NURBS surface. From the relationship between ru,rv,rw, it was known that rw is perpendicular to the plane formed by ru and rv, so the sum of ru⋅rw and rv⋅rw was 0. Then, we could obtain:

{S⋅ru=Bru2+Crv⋅ruS⋅rv=Crv2+Brv⋅ru(2)

The equation was solved as:

B=S⋅ru−Brv⋅ruru2C=S⋅rvru2−S⋅ru(rv⋅ru)ru2rv2−(ru⋅rv)2(3)

After calculation, the node vector was updated for calculation by:

{u=u−Bv=v−C(4)

Finally, the calculation error ε was given as:

|S⋅ru|≤ε,|S⋅rv|≤ε(5)

When the given conditions were met, the point P1 is the desired vertical foot point, and the iteration process was terminated. Otherwise, the iteration continued until the conditions were met. A method [50,51] was used to calculate the distance function SN. It was supposed that there were two NURBS bodies A and B in contact, and two boundary points C(1) and C(2) that belonged to the two bodies respectively. Cvp(2) is the nearest projection point of the point C(1) on the contact boundary of object B. As C(2)(vp1,vp2) was set as the NURBS surface of the contact boundary of object B, then SN was be expressed as:

SN=C(1)−C(2)(vp1,vp2)(6)

From the given conditions it was obtained that:

{f1=[C(1)−C(vp1,vp2)]⋅C1′(vp1,vp2)=0f2=[C(1)−C(vp1,vp2)]⋅C2′(vp1,vp2)=0(7)

Newton’s method was used to solve Eq. (7), and the corresponding iterative equation was:

H(k)[Δvp1Δvp2](k+1)=−[[C(1)−C(vp1,vp2)]⋅C1′(vp1,vp2)[C(1)−C(vp1,vp2)]⋅C2′(vp1,vp2)](8)

where, k is the number of iterations and H is the Jacobian matrix:

H=[∂(f1)∂vp1∂(f1)∂vp2∂(f2)∂vp1∂(f2)∂vp2](9)

3.2 Treatment of Constraints and Loads

When dealing with contact constraints, numerical solutions should be carried out. The widely used penalty function method is used to calculate the contact force [52], however, this method does not increase the degree of freedom in the numerical solution and does not change the scale of the matrix [53]. The penalty factor introduced by the penalty function was multiplied by the penetration of the contact interface to calculate the contact force. In fact, the contact problem was transformed into a material problem to solve, and the specific calculation formula was expressed as:

PN=εs(10)

where, ε is the penalty factor and s is the amount of invasion. The constraint processing was to transfer the constraint conditions into a functional solution. The specific functions were expressed as:

∏=∏e+∏c(11)

∏c=12∫ΓcεSN2dΓc(12)

The equation was rewritten by variation:

δ∏c=∫ΓcεδSNSNdΓc(13)

According to Eq. (6), Eq. (13) was rewritten as:

δ∏c=∫ΓcC(1)−C(2)(vp1,vp2)(δC(1)−C(2)(vp1,vp2))dΓc(14)

Since the contact surfaces of A and B were composed of NURBS surfaces, we could set:

C(1)=∑I=1NP1NI(1)qI(1),C(2)=∑J=1NP2NJ(2)qJ(2)(15)

where, N(i) is a shape function, q(i) are the coordinates of NURBS control points, and NPi is the number of control points on the contact boundary of two objects. To their variations:

δC(1)=∑I=1NP1NI(1)δqI(1),δC(2)=∑J=1NP2NJ(2)δqJ(2)(16)

Assuming that the set of Gauss points on the contact area was defined as G(A), the i-th Gaussian point of it is C(1)(ui). For one Gaussian point on A, there would be one corresponding nearest projection point on B, so Eq. (14) could be rewritten as:

δ∏c=∑i∈G(1)εωiJ(1)(ui)[C(1)−C(2)(vp1,vp2)][∑I=1NP1NI(1)qI(1)−∑J=1NP2NJ(2)qJ(2)](17)

where, ωi is the integral weight, J(1)(ui) is the Jacobian matrix of the isoparametric transformation, J(1)(ui)=‖C1(1)×C2(1)‖. Therefore, Eq. (17) could be rewritten as:

δ∏c=−∑I=1NP1RI(1)δqI(1)−∑J=1NP2RJ(2)δqJ(2)(18)

In this way, the projection of the contact force at the control point was obtained as:

{RI(1)=∑i∈G(1)εωiJ(1)(ui)NI(1)[C(2)(vp1,vp2)−C(1)]RJ(2)=∑i∈G(1)εωiJ(1)(ui)NJ(2)[C(1)−C(2)(vp1,vp2)](19)

The three-dimensional contact stiffness matrix was expressed as:

KJC=[KJC(11)KJC(12)KJC(21)KJC(22)](20)

where,

KJCij=∂RI(i)∂qJ(j)=∂RI(i)∂qJ(j)|vp+∂R(i)∂vp1⋅∂vp1∂qJ(j)+∂R(i)∂vp2⋅∂vp2∂qJ(j)(21)

Suppose that:

t1(2)=C,vp1(2)=∂C(2)(vp1,vp2)∂vp1t2(2)=C,vp2(2)=∂C(2)(vp1,vp2)∂vp2(22)

t1(2),t2(2) are the tangent vectors of the control points.

{f1=[C(1)−C(2)(vp1,vp2)]⋅t1(2)=0f2=[C(1)−C(2)(vp1,vp2)]⋅t2(2)=0(23)

Newton’s method was used to solve the differential equation, and the iterative formula was:

H(k)[Δvp1Δvp2](k+1)=−[[C(1)−C(vp1,vp2)]⋅t1(2)[C(1)−C(vp1,vp2)]⋅t2(2)](k)(24)

where H is the Jacobian matrix:

H=[∂(f1)∂vp1∂(f1)∂vp2∂(f2)∂vp1∂(f2)∂vp2]=[−t1(2)⋅t1(2)+S⋅C11(2)(vp1,vp2)−t1(2)⋅t2(2)+S⋅C12(2)(vp1,vp2)−t2(2)⋅t1(2)+S⋅C21(2)(vp1,vp2)−t2(2)⋅t2(2)+S⋅C22(2)(vp1,vp2)](25)

The variation of Eq. (23) was obtained as:

{[C(1)−C(2)(vp1,vp2)]⋅δt1(2)+[δC(1)−δC(2)(vp1,vp2)]⋅t1(2)=0[C(1)−C(2)(vp1,vp2)]⋅δt2(2)+[δC(1)−δC(2)(vp1,vp2)]⋅t2(2)=0(26)

To obtain the stiffness matrix [54], it was necessary to further deduce and perform first-order Taylor expansion on vp1,vp2 in the following three items:

{C(2)=C(2)(vp1,vp2)+∂C(2)(vp1,vp2)∂vp1|vp1(v−vp2)+∂C(2)(vp1,vp2)∂vp2|vp2(v−vp2)=C(2)(vp1,vp2)+t1(2)(v−vp1)+t2(2)(v−vp2)C,vp1(2)=C,vp1(2)(vp1,vp2)+∂C,vp1(2)(vp1,vp2)∂vp1|vp1(v−vp2)+∂C,vp1(2)(vp1,vp2)∂vp2|vp2(v−vp2)=t1(2)+C,11(2)(vp1,vp2)(v−vp1)+C,12(2)(vp1,vp2)(v−vp2)C,vp2(2)=C,vp2(2)(vp1,vp2)+∂C,vp2(2)(vp1,vp2)∂vp1|vp1(v−vp2)+∂C,vp2(2)(vp1,vp2)∂vp2|vp2(v−vp2)=t2(2)+C,21(2)(vp1,vp2)(v−vp1)+C,22(2)(vp1,vp2)(v−vp2)(27)

Eq. (27) was obtained by variation as:

δC(2)=δC(2)(vp1,vp2)−t1(2)δvp1−t2(2)δvp2δC,1(2)=δt1(2)−C,11(2)(vp1,vp2)δvp1−C,12(2)(vp1,vp2)δvp2δC,2(2)=δt2(2)−C,21(2)(vp1,vp2)δvp1−C,22(2)(vp1,vp2)δvp2(28)

Substituting Eq. (28) into Eq. (26), we obtained:

{[−t1(2)⋅t1(2)+S⋅C,11(2)(vp1,vp2)]δvp1+[−t2(2)⋅t1(2)+SC,12(2)(vp1,vp2)]δvp2=SδC,1(2)+[−δC(1)+δC(2)]t1(2)[−t1(2)t2(2)+SC,21(2)(vp1,vp2)]δvp1+[−t2(2)t2(2)+SC,22(2)(vp1,vp2)]δvp2=SδC,2(2)+[−δC(1)+δC(2)]t2(2)(29)

Combined with the Jacobian matrix H, we obtained:

H[δvp1δvp2]=−[S⋅δC,1(2)+[δC(1)−δC(2)]⋅t1(2)S⋅δC,2(2)+[δC(1)−δC(2)]⋅t2(2)](30)

The NURBS discretization of the above formula was obtained:

[δvp1δvp2]=−[H−1](S[∑J=1NP2NJ,1(2)δqJ(2)∑J=1NP2NJ,2(2)δqJ(2)]+[∑I=1NP1NI(1)δqI(1)−∑J=1NP2NJ(2)δqJ(2)][t1(2)t2(2)])(31)

So, the terms in Eq. (21) were gotten as:

{[∂vp1∂qI(1)∂vp2∂qI(1)]=−[H](NI(1)[t1(2)t1(2)])[∂vp1∂qJ(j)∂vp2∂qJ(j)]=[H](NJ(2)[t1(2)t1(2)]−[NJ,1(2)⋅SNJ,2(2)⋅S])(32)

{∂RI(1)∂qI(1)=−εωiJ(1)(ui)NI(1)NI(1)∂RI(1)∂qJ(2)=εωiJ(1)(ui)NI(1)NJ(2)∂RJ(2)∂qI(1)=εωiJ(1)(ui)NJ(2)NI(1)∂RJ(2)∂qJ(2)=−εωiJ(1)(ui)NJ(2)NJ(2)(33)

{[∂RI(1)∂vp1∂RI(1)∂vp2]=εωiJ(1)(ui)NI(1)[t1(2)t2(2)][∂RJ(2)∂vp1∂RJ(2)∂vp2]=−εωiJ(1)(ui)NJ(1)[t1(2)t2(2)](34)

Then, the contact matrix KJC was obtained:

KJC(11)=∑i∈G(1)GiNI(1)NI(1)[A(2)−I]KJC(12)=∑i∈G(1)Gi(NI(1)NJ(2)[I+A(2)]−NJ(2)B(2))KJC(21)=∑i∈G(1)GiNI(1)NJ(2)[A(2)+I]KJC(22)=∑i∈G(1)Gi(NI(2)NJ(2)[−I−A(2)]+NJ(2)B(2))(35)

where,

Gi=εωiJ(1)A(2)=[t1(2)t2(2)][H−1][t1(2)t2(2)]B(2)=[t1(2)t2(2)][H−1][NJ,1(2)⋅SNJ,2(2)⋅S](36)

After obtaining the stiffness matrix, the equilibrium equations of the contact system were obtained, which can be expressed in matrix form as:

[K1+KJC(11)KJC(12)KJC(21)K2+KJC(22)][u(1)u(2)]=[f(1)f(2)](37)

where, u(1) and u(2) represent the displacement vectors of the control points on object A and object B, f(1) and f(2) represent the external force loads applied on A and B, respectively. The calculation of K1,K2 in the equation was similar to that of FEA, which could be obtained by using the principle of minimum potential energy.

Kq=F(38)

Then, the element stiffness matrix was written as:

Ke=∫VeBTDBdV(39)

After the element stiffness matrix was obtained, the overall stiffness matrix K=∑Ke was assembled according to the sequence of control points, where, D is the elastic matrix and B is the geometric matrix.

D=E(1−μ)(1+μ)(1−2μ)⋅[1μ1−μμ1−μ000μ1−μ1μ1−μ000μ1−μμ1−μ10000001−2μ2(1−μ)0000001−2μ2(1−μ)0000001−2μ2(1−μ)](40)

where, E is the elastic modulus and μ is Poisson’s ratio.

B=[dN1dxdN2dx…dNidxdN1dydN2dy…dNidydN1dzdN2dz…dNidzdN1dydN2dy…dNidydN1dxdN2dx…dNidxdN1dydN2dy…dNidydN1dzdN2dz…dNidzdN1dzdN2dz…dNidzdN1dxdN2dx…dNidx](41)

where, Ni is the NURBS basis function, i is the total number of control points on a unit in a NURBS volume.

Although isoparametric elements are used in isoparametric geometry, the mapping relationships among parent elements coordinate, physical coordinate, and parameter coordinate need to be carried out in IGA [55]. The parent element refers to the Gaussian integral domain. As shown in Fig. 3, ϕ refers to the mapping between the physical unit coordinate and parameter space coordinate, and ϕ1 refers to the mapping between the parameter coordinate and parent unit coordinate.

Figure 3: Mapping among domains

The determinant of Jacobian transformation corresponding to the transformation from parameter coordinate to parent unit coordinate is:

|J1|=(ζi+1−ζi)(ηj+1−ηj)(ξi+1−ξi)/8(42)

The Jacobian transformation matrix corresponding to the transformation from physical unit coordinate to parameter space coordinate is:

J2=[dxdζdxdηdxdξdydζdydηdydξdzdζdzdηdzdξ](43)

Therefore, the Jacobian determinant of the mapping from the parent unit to the physical unit was obtained as:

|J|=|J1||J2|(44)

The expression of the element stiffness matrix was:

Ke=∑r=1nr∑k=1mk∑j=1mj∑i=1miωijkBTDB|J|(45)

where, nr is the total number of units of the three-dimensional model. mi,mj,mk are the number of Gaussian integral points in the direction of i,j,k. ωijk is the weight of Gaussian point. Finally, Eq. (39) was solved by Gaussian integral.

3.3 Hertz Contact Example Validation

For the Hertz contact problem, this subsection gives the geometric settings and boundary conditions of the Hertz contact model, calculates the analytical expression of the Hertz contact problem, and compares it with the FEA models with different mesh densities to verify the correctness of the equivalent geometric contact algorithm.

3.3.1 Theoretical Value of the Hertz Contact Model

In 1881, Hertz obtained the formulas required for calculating the contact problem according to the mathematical linear elasticity method and gave the following assumptions: Each object is regarded as an elastic, uniform and isotropic half space, which is loaded on a small elliptical region of its smooth surface. Assuming that the surface has no friction, its contact area must be smaller than the surface radius and object size, and the geometry of the contact interface is described by quadratic polynomial. The Hertz contact model is shown in Fig. 4.

Figure 4: Hertz contact model

The analytical solution of the half width b of the contact area was obtained from Eq. (46).

b=4Fπl(1−μ12)/E1+(1−μ22)/E2(1/r1)+(1/r2)(46)

As evident from Eq. (46), the value of half width b depends on the applied concentrated force F, the Poisson’s ratio of the upper and lower cylinders are μ1 and μ2, and the corresponding elastic modulus are E1 and E2, and the radius of the cylinder are r1 and r2, and its length in the x-direction is L. The stress distribution of half-width b in the contact area was calculated from Eq. (47).

p(y)=2Fπlb2b2−y2=12r1+12r21−μ12E1+1−μ22E22F1−μ12E1+1−μ22E2πl(12r1+12r2)−y2(47)

Considering the symmetry of the model, a quarter model could be established to save on calculation costs. The radius of the two cylinders was 2 mm, the thickness was 5 mm, and the Neumann boundary condition F was 1000 N. To reduce the influence of simplifying the mechanical model, the constant Dirichlet boundary condition v = 0.001125 was used to replace the Neumann boundary condition F at the top of the cylinder [49]. The penalty function method was used for contact analysis, and the penalty value is 1000E. According to Eqs. (46) and (47), the Hertz contact theoretical value was calculated as: the contact half width b = 0.0237 mm and the maximum contact stress value was 735 MPa. The Hertz contact analysis model and calculation parameters are shown in Fig. 5.

Figure 5: Hertz contact analysis model

3.3.2 Hertzian Contact Values Calculated by IGA and FEA

To fully verify the accuracy of the isogeometric contact analysis, the isogeometric Hertz contact analysis and FEA were carried out concurrently. The advantages of the concurrent analysis were mainly reflected in that it allowed a visual determination of whether the 3D isogeometric contact theory and procedure could pass the Hertz contact test, and additionally, by adjusting the computational unit densities and comparing it with FEA, the differences or advantages between the results of IGA and FEA under different computational unit densities were explored. The isogeometric Hertz contact model used NURBS body with order 2 in u, v and w directions. To consider the influence of calculation element density on calculation accuracy, five NURBS body models with different computational unit densities were obtained by adjusting the number of refinements, and the node vectors in the u−direction (or v−direction or w−direction) of the five NURBS body models were {0 0 0 1 1 1}, {0 0 0 0.5 0.5 1 1 1}, {0 0 0 0.25 0.5 0.5 0.75 1 1 1 }, {0 0 0 0.125 0.25 0.375 0.5 0.5 0.625 0.75 0.875 1 1 1}, {0 0 0 0.0625 0.125 0.1875 0.25 0.3125 0.375 0.4375 0.5 0.5 0.5625 0.625 0.6875 0.75 0.8125 0.875 0.9375 1 1 1}. Due to the tensor product property, it was difficult to achieve local subdivision of NURBS volume. To ensure unity, the global refinement modeling method was also used for the FEA as a comparison. The mesh sizes of the five finite element models were 0.25, 0.2, 0.15, 0.1, and 0.05 mm, respectively.

The selection of FEA contact settings and solvers had a significant impact on the analysis results. The commercial finite element software selected in this article was Ansys Workbench. For contact settings, the contact surfaces of the upper and lower half cylinders were selected for contact and target surfaces, respectively; select Contact Type as Frictionless. The penetration tolerance was set to 0.001, the penalty function value was 10E, and the contact offset value was 0 mm. The “Contact Geometry Correction” and “Target Geometry Correction” options were disabled. The rest were default settings. At the Solver setup, the Iteration method selection was the Newton-Raphson method; the residual value was 10E-6, and the rest were default settings.

Generally, the more accurate calculation results were due to the large number of calculation units. Therefore, the IGA results with the node vector {0 0 0 0 0.0625 0.125 0.1875 0.25 0.3125 0.375 0.4375 0.5 0.5 0.5625 0.625 0.6875 0.75 0.8125 0.875 0.9375 1 1 1} were compared with the FEA results with the mesh size of 0.05 mm. The results are shown in Fig. 6, where Fig. 6a shows the results of equivalent stress distribution, Fig. 6b shows the results of stress distribution in Y–direction, and Fig. 6c shows the results of contact stress distribution.

Figure 6: Comparison of IGA and FEA Hertz contact model analysis results

According to the Hertz contact theory value, the maximum contact stress of the two columns was 735 MPa. From the analysis results in Fig. 6c, the maximum contact stress calculated by FEA was 704.88 MPa, and the maximum contact stress calculated by IGA was 719.27 MPa. The results of IGA and FEA were close to the theoretical value in the case of dense calculation elements. Figs. 6a and 6b show that the maximum equivalent force and maximum Y-directional stress in the FEA results were 475.57 MPa and −857.03 MPa, while the maximum equivalent force and maximum Y-directional stress in the IGA results were 482.43 and −866.32 MPa, respectively. In terms of stress distribution, the results of IGA and FEA had approximately the same trend of stress distribution.

To compare the effect of computational unit densities on the analysis results, the contact stresses of isogeometric and finite element models with different mesh densities were studied, the boundary effect on the contact stress was also eliminated, and the contact stress at the contact position of the cross-section in the half thickness direction was selected and compared with the Hertz contact theory value. Fig. 7 shows the maximum contact stress values calculated by IGA and FEA at different computational unit densities and compares the maximum contact stress value with the maximum Hertz contact theoretical value.

Figure 7: Convergence of maximum contact pressure between FEA and IGA

As shown in Fig. 7, with the increase in calculation elements, the IGA and FEA results converged to the theoretical value. Compared with the FEA, the isogeometric contact analysis program obtained higher calculation accuracy with fewer calculation elements, and the increase of calculation elements seemed to have less impact on the results of IGA. The isogeometric models with five computational unit densities were compared with the element and node data of the finite element model, as presented in Table 1.

Under the same calculation accuracy, the number of elements used in the isogeometric contact analysis was less. Taking the 4× refinements as an example, the number of IGA units was only 34,992, while the number of FEA units was 290,433, which was eight times the number of IGA units. The accuracy of IGA was less affected by the density of computational units and converged more easily than the finite element method. Fig. 8 shows the comparison between the results of the geometric contact stress distribution of the refined 2, 3, and 4 contact width areas and the Hertz contact theoretical value, which also showed that the accuracy of the geometric analysis was less affected by the element density.

Figure 8: Contact stress distribution in contact width region

In summary, the IGA contact algorithm in this paper passed the Hertz contact test. Compared with the FEA method, under global refinement, the IGA contact algorithm had fewer nodal degrees of freedom and computational units with the same computational accuracy, and the analytical accuracy was less affected by the computational unit density, which was easier to converge.

4  Isogeometric Contact Analysis for Gears

Compared with the traditional FEA method, the IGA method has the greatest advantage of using the same geometric model representation in modeling and analysis, which can eliminate the discrete error and effectively avoid the discrete approximation problem of the finite element model. Additionally, the boundary of the IGA model can be at least C(1) continuous, which can theoretically yield more accurate calculation results. At present, the models of IGA mostly use NURBS splines. The mathematical expression of a NURBS body is similar to that of curves and surfaces, and is expressed by multiplying the primary functions in three directions, as shown in Eq. (48). Assuming that the degree of the basis function in the three directions of u, v, and w is p, q, and r, the NURBS body can be defined as:

{V(u,v,w)=∑i=1n∑j=1m∑k=1lRi,j,k(u,v,w)Pi,j,kRi,j,k(u,v,w)=Ni,p(u)Nj,q(v)Nk,r(w)ωi,j,k∑i=1n∑j=1m∑k=1lNη,p(u)Nς,q(v)Nζ,r(w)ωη,ς,ζ(48)