## 1 Introduction  
The gearbox is the key unit of modern machinery and has been widely used in all kinds of engineering situations. Due to the complex structure and tough running environments of large-scale machinery, key gearbox components emerge easily to local defects, which brings about fatal accidents and generates major disruptions [1; 2]. Therefore, timely detection of faults residing in rolling bearings and gears are valuable for ensuring the safe and stable operation of mechanical equipment [3; 4; 5; 6]. Meshing stiffness is the ability of the teeth to resist deformation during meshing, and it is the key factor of gear vibration response. Therefore, it is of great significance to study the time-varying meshing stiffness. In order to better study the various working conditions of rotating equipment, it is necessary to study its dynamic modeling. Scholars have also done a lot of research in the field of dynamics simulation. When calculating the meshing stiffness, the analytical method is widely used due to its efficiency and accuracy when comparing with FEM (finite element method). [PERSON] [7] studied the meshing stiffness when accounting for loading and stressing distribution on helical gears; the results matched well with the Finite Element Analysis (FEA)method. [PERSON] [8] used the FEA to calculate the torsional mesh stiffness of gears in mesh when there is a crack in the tooth root. [PERSON] [9] used a 2D FEM model to calculate the meshing stiffness in different crack depths, In this model, the actual meshing force was replaced by the uniformly distributed load. By drawing the boundary of the root stress, the effect of tooth thickness could be calculated, and then the meshing stiffness could be derived. The accuracy of the 2D FEM method was proven by the photoelastic experiment method by [PERSON] and [PERSON] [10]. The 3D FEM method was established by [PERSON] [11]. By putting a uniformly distributed load in the tooth profile, the displacement in the meshing direction was extracted, and the meshing stiffness could then be derived by calculating the ratio between the meshing force and the displacement. Above all, the 2D FEM method is a widely used method, matches well with the experimental result, and could be the right reference for the analytical method. [PERSON] et al. [12] analyzed the dynamic model signal by their method, and the results were consistent with the experimental signals, which proves the method is effective. Analytical methods are widely used methods due to the high efficiency when compared with FEM methods. The Silurian method was proposed by [PERSON] [13]. In this method, the gear profile is simplified by regarding the gear as a complex of a trapezoid and a rectangular, so it is far away from the actual gear profile. The relationship between meshing stiffness and potential energy was derived by [PERSON] and [PERSON] [14]. On the basis of [PERSON] and Lin, a time-varying meshing stiffness calculation method that accounts for the Hertzian, bending, shear and axial compressive energies was proposed by [PERSON] [15], but the fill-foundation stiffness has not been considered. The theoretically precise calculation method of fillet foundation deflection, which was verified by FEM, was originally proposed by [PERSON] [16]. Then, the effect of different types of cracks on the dynamic responses was investigated by [PERSON] [17] by taking the fillet-foundation stiffness into account. On this basis, [PERSON] et al. [18] proposed a quantitative trend fault diagnosis method of a coupling gearbox based on Sparsogram and Lempel-Ziv. The performance is much better on the quantitative diagnosis of cracked components. Nevertheless, the tooth profile between the base circle and the root circle had rarely been considered. To solve this problem, an improved mesh stiffness calculating method was proposed by [PERSON] [19] by considering the misaligned behavior between the base circle and the gear root circle; however, the expression of equations should be changed by the actual tooth numbers.