Gear pitting kinetic model solving method based on microscopic surface topography

Abstract

The invention belongs to the technical field of gear transmission system fault diagnosis, and particularly relates to a gear pitting kinetic model solving method based on microscopic surface topography. Comprising the following steps: establishing a tooth surface morphology and pitting defect model; the tooth surface morphology of the selected gear is sampled, and nonlinear excitation parameters in a gear transmission system are collected; a gear tooth surface morphology mathematical model based on the fractal dimension D and the fractal roughness G is established, and a mixed elastohydrodynamic lubrication coupling analysis mathematical model of the pitted gear under the rough surface is established based on the relation that the curvature radius, the dynamic load, the oil film entrainment speed and the like change along with the change of the meshing position under different surface roughness. And solving the non-linear dynamic mathematical model of the gear with the rough surface. According to the method, the influence of the surface roughness on the rigidity attenuation and pitting extension of the gear can be accurately quantified, the nonlinear dynamic behavior of the gear system under the complex lubrication condition is disclosed, and an analysis method support is provided for the pitting resistance design and life prediction of the gear.

Description

Technical Field

[0001] The present invention belongs to the technical field of gear transmission system fault diagnosis, and in particular relates to a method for solving a gear pitting dynamics model based on microscopic surface morphology. Background Art

[0002] During meshing transmission, gear tooth surface roughness is closely linked to tooth backlash and time-varying mesh stiffness. The height of the asperities on a rough tooth surface affects the size of the backlash during meshing, which in turn affects the normal contact stiffness of the meshing interface, ultimately affecting the dynamic behavior of the gear system. Tooth surface roughness also affects gear lubrication. During gear operation, lubricant acts as a hydrodynamic wedge, forming an oil film of a certain thickness between the meshing teeth. Due to the characteristics of the gear, the load and other operating conditions acting on the oil film vary, causing deformation and changes in film stiffness. Furthermore, the viscous nature of the lubricant causes film damping, leading to energy dissipation in the gear system during operation. This complicates the motion of the gear system under lubricated conditions. Tooth surface roughness significantly affects the oil film thickness, affecting gear film damping, film stiffness, and friction coefficient, ultimately impacting the nonlinear dynamic characteristics of the gear system. As gears continue to wear in the stable wear phase over time, their surface roughness gradually increases, significantly affecting the lubrication state and dynamic response of the gear system. Traditional dynamic models often assume smooth tooth surfaces and ignore the impact of roughness on stiffness degradation and pitting expansion, resulting in insufficient analytical accuracy. Summary of the Invention

[0003] The purpose of the present invention is to provide a method for solving the gear pitting dynamics model based on microscopic surface morphology. The method performs three-dimensional morphology characterization of the tooth surface roughness based on fractal theory, establishes a tooth surface contact stiffness correction model considering micropitting defects, and combines the hybrid elastohydrodynamic lubrication characteristics to analyze the influence of oil film stiffness, oil film damping and time-varying friction coefficient on the dynamic response of the system. The present invention can accurately quantify the influence of surface roughness on gear stiffness attenuation and pitting expansion, reveal the nonlinear dynamic behavior of the gear system under complex lubrication conditions, and provide analytical method support for gear anti-pitting design and life prediction.

[0004] To achieve the above objectives, the present invention adopts the following technical solutions.

[0005] A method for solving a gear pitting dynamics model based on microscopic surface morphology includes steps 1 to 4. The specific contents of each step include:

[0006] Step 1. Establish the tooth surface morphology and pitting defect model; specifically divided into the following sub-steps:

[0007] 1a. Sample the tooth surface topography of the selected gear and collect the non-linear excitation parameters inside the gear transmission system;

[0008] 1b. Establish a mathematical model of the gear tooth surface topography based on the fractal dimension D and the fractal roughness G, which can be expressed as: ;

[0009] where, is the power spectrum, ω is the spatial frequency, and γ is the constant related to the spatial frequency;

[0010] The calculation methods of the fractal dimension D and the fractal roughness G are as follows: Plot the logarithm of the spatial frequency lgω on the abscissa and the logarithm of the power spectral density lgP(ω) on the ordinate to obtain the corresponding slope and intercept, and calculate the corresponding fractal dimension D and fractal roughness G;

[0011] 1c. Establish a mathematical model of gear tooth surface deformation

[0012] Use the W-M function combined with the fractal roughness to express the gear tooth surface deformation, and the expression of the fractal rough surface is obtained as: ;

[0013] In the formula, z i (t), where i = p or g, respectively refer to the fractal rough surfaces of the driving wheel and the driven wheel; λ is the scale coefficient, D is the fractal dimension, 1 < D < 2, n is the frequency index, is the maximum frequency index; and respectively represent the fractal rough surfaces of the driving wheel and the driven wheel, represents the clearance between the two tooth profiles without considering the surface topography, is the scale coefficient, and t is the time;

[0014] L g is the sampling length of the fractal rough surface, expressed as ;

[0015] where, R m =R p R g / (R p +R g ) is the comprehensive curvature radius of the meshing pair, R p is the curvature radius of the driving wheel, R g is the curvature radius of the driven wheel, F N is the meshing force; E is the elastic modulus;

[0016] Furthermore, the unilateral tooth side clearance of the gear considering the surface topography is obtained as:

[0017]

[0018] Where, is the single-sided tooth backlash of the j-th gear pair, is the initial single-sided tooth backlash without considering the surface topography; It refers to the fractal dimension of the driving wheel and the driven wheel. It refers to the fractal dimension of the driving wheel and the driven wheel;

[0019] use Represents the projection length of the relative displacement of the jth pair of gears on the meshing line, then the deformation between the gear pairs is expressed by the piecewise function Expressed as:

[0020]

[0021] Meshing force between gear pairs j Expressed as: ;in represents the meshing stiffness of the j-th gear pair; is the derivative of the projected length of the relative displacement of the jth gear pair on the meshing line;

[0022] Step 2: Establish a time-varying mesh stiffness correction model

[0023] 2a. Establish a mathematical solution model for contact stiffness optimization considering surface topography

[0024] Consider the gear teeth as a nonlinear cantilever beam with variable cross-section. Based on step 1 and fractal theory, the stiffness of the contact interface is considered as the combination of the stiffness of all asperities. Based on the contact area and critical length, the contact stiffness considering the surface topography can be obtained by using quadratic integration. :

[0025] ;

[0026] Where D and G are the fractal dimension and fractal roughness of meshing, respectively, and ; ; represents the stiffness of a single asperity and ; is the actual contact area of ​​a single asperity; is the critical length of the asperity and , is Poisson's ratio, is the yield strength; is the distribution function of the asperity, Represents the maximum contact area of ​​a single micro-convex body, La is the diameter of the bottom surface of the micro-convex body, , represents the surface coefficient; D p is the fractal dimension of the driving wheel in the gear pair, D g is the fractal dimension of the driven wheel, Gp With G g are the fractal roughness of the driving wheel and the driven wheel in the gear pair respectively;

[0027] 2b. Establish a mathematical model for optimizing meshing stiffness considering surface topography, including:

[0028] (2B1) Mathematical model when the base circle diameter is smaller than the root circle radius

[0029] In the mathematical model when the base circle diameter is smaller than the root circle radius, define F a With F b The meshing force F Nj The components of force in the horizontal and vertical directions; d(y) represents the distance between the meshing point and the base circle in the direction of tooth height; h x represents the distance from the base circle x to the center line of the gear tooth; h(y) represents the distance between the meshing point and the center line of the gear tooth in the direction of tooth thickness; h f Indicates the distance from the intersection of the root circle and the tooth profile to the center line of the gear tooth; h b Indicates the distance from the starting point of the involute on the base circle to the center line of the gear tooth; Indicates the meshing force F Nj With F b Angle; Indicates the angle between the starting point of the involute on the base circle and the center line and the center line of the gear tooth; Indicates the angle between the starting point of the involute on the root circle and the center line of the gear tooth; Indicates the angle between the tooth centerline and the line connecting the base circle tangent points; Indicates the pitch circle pressure angle, represents the root circle pressure angle; r b Indicates the radius of the gear base circle; r f Indicates the tooth root circle radius;

[0030] Then the bending potential energy Ub can be expressed as: ;

[0031] The shear potential energy Us can be expressed as: ;

[0032] The radial compression potential energy Ua can be expressed as: ;

[0033] Where, E, Represent the gear elastic modulus and shear modulus, I x A represents the moment of inertia of the tooth section at a distance x from the base circle in the tooth height direction. x represents the cross-sectional area of ​​the above section;

[0034] The parameter derivation formula is as follows:

[0035] ;

[0036] ;

[0037] ;

[0038] ;

[0039] ;

[0040] ;

[0041] Substituting the parameter derivation formula into the potential energy expression, we can obtain the Hertz contact stiffness , bending stiffness , shear stiffness The relationship with parameters:

[0042] ;

[0043] ;

[0044] ;

[0045] ;

[0046] ;

[0047] ;

[0048] ;

[0049] (2B2) Mathematical model when the base circle radius is greater than the root circle radius

[0050] When the base circle radius is larger than the root circle radius, the gear tooth deformation between the root circle and the base circle is taken into account in the calculation model, and the result is:

[0051] ;

[0052] ;

[0053] ;

[0054] ;

[0055] ;

[0056] ;

[0057] in Refers to the distance between the pitch circle and the tooth root position. It refers to the distance from a point on the transition curve from tooth root to involute to the center line of the gear tooth. Indicates the tooth root transition curve radius;

[0058] Further we get:

[0059] ;

[0060] in ;

[0061] ;

[0062] ;

[0063] (2B3) Mathematical model of internal gear

[0064] In the mathematical model of internal gears, define D r is the root of the inner gear ring, F is the meshing point, X is D r Any point on the line connecting points X and F, x represents the distance between points X and F in the direction of the center line, h represents the distance between points X and F in the direction of the center line, x is the distance from point X to the center line of the inner ring gear teeth, d is the distance between the meshing point and the tooth root in the center line direction, It is the angle between the tangent line of the base circle through point X and the center line of the tooth groove. The angle between the tangent line of the base circle at point F and the center line of the tooth groove is, for The angle between the component perpendicular to the center line and the component perpendicular to the center line is calculated by transforming the corresponding calculation method according to the angle relationship in the cantilever beam:

[0065] ;

[0066] ;

[0067] ;

[0068] in ,

[0069] ;

[0070] ;

[0071] ;

[0072] ;

[0073] It refers to the base circle radius; refers to the tooth root radius; B refers to the deformation coefficient; It refers to the angle between the starting point of the involute on the root circle and the center line of the gear tooth;

[0074] Step 3: Establish a mathematical model for the coupled analysis of hybrid elastohydrodynamic lubrication of pitting gears under rough surfaces

[0075] Based on the relationship between the curvature radius, dynamic load, oil film entrainment speed and other factors as the meshing position changes under different surface roughness, a mathematical model for the hybrid elastohydrodynamic lubrication coupling analysis of pitting gears under rough surfaces is established, including:

[0076] (3a) The lubrication model of the gear system under stable working condition is expressed as:

[0077] ;

[0078] In the lubrication model of the gear system under stable working conditions, x represents the coordinate in the same direction as the main speed, ρ represents the density of the lubricating oil, h represents the thickness of the lubricating oil film, η represents the viscosity of the lubricating oil, p represents the oil film pressure on the lubricating oil, and u o1 、u o2 Represent the velocities of the upper and lower surfaces along the x direction, u o represents the entrainment velocity. The physical meaning of the lubrication model under the stable working state of the gear system is the balance of flow. The left side of the model is the pressure flow formed by the pressure gradient of the lubricating oil film along the x direction, and the rightmost side of the model is the shear flow caused by the surface velocity.

[0079] (3b) The mathematical model of film thickness in the mixed elastohydrodynamic lubrication process is expressed as:

[0080] ;

[0081] In the film thickness mathematical model of the mixed elastohydrodynamic lubrication process, h(x) is the lubricating oil film thickness at x, h0 is the central film thickness, R is the instantaneous curvature radius, υ(x) is the elastic deformation at x, and s r (x) is the surface roughness at position x;

[0082] (3b) Mathematical model of viscosity and density in mixed elastohydrodynamic lubrication process

[0083] The viscosity mathematical model of the mixed elastohydrodynamic lubrication process is obtained by fitting the experimental results and is expressed as:

[0084] ;

[0085] In the viscosity mathematical model of the mixed elastohydrodynamic lubrication process, T0 is the initial temperature, is the lubricant viscosity at T=T0, is the pressure viscosity coefficient, is a constant;

[0086] The density mathematical model of the mixed elastohydrodynamic lubrication process is expressed as:

[0087] ;

[0088] In the density mathematical model of the mixed elastohydrodynamic lubrication process, is the lubricant density when p=0, T=T0, is the density-temperature coefficient;

[0089] Step 4: Solve the nonlinear dynamic mathematical model of gears with rough surfaces

[0090] Based on steps 1 to 3, the solution process for the dynamic characteristics of the gear transmission with pitting corrosion considering the rough tooth surface includes:

[0091] (4a) Use an optical profilometer to detect gear tooth roughness and confirm the basic parameters and nonlinear excitation parameters of the gear;

[0092] (4b) Substituting the calculated dynamic internal excitation into the derived corresponding mathematical model for iterative solution to solve the relative meshing displacement and relative meshing speed of the system. At the beginning of the iteration, the friction torque is set to 0;

[0093] (4c) After the relative meshing displacement and relative meshing speed are calculated, the relative entrainment speed, meshing force, and equivalent curvature radius are determined based on the basic parameters and nonlinear excitation parameters of the gear, the relative meshing displacement, and the relative meshing speed. The kinematic and load parameters are substituted into the hybrid elastohydrodynamic lubrication model every Δt to iterate and obtain the hybrid elastohydrodynamic lubrication characteristics at that moment. The hybrid elastohydrodynamic lubrication characteristic data within one meshing cycle are calculated.

[0094] (4d) Substitute the obtained time-varying hybrid elastohydrodynamic lubrication characteristic data into the dynamic mathematical model for the next iteration until the model converges. The iterative convergence condition is:

[0095] ;

[0096] Where Q is the number of iterations, and γ is taken as all meshing positions calculated within one meshing cycle;

[0097] is the gear meshing force, It refers to time, Refers to the gear that bears the force. Refers to the number of gears.

[0098] A further improvement or preferred implementation of the aforementioned method for solving the gear pitting dynamics model based on microscopic surface morphology, wherein the nonlinear excitation parameters in step (1a) include three-dimensional surface morphology, two-dimensional depth morphology, profile (depth, width, curvature, angle), and surface roughness.

[0099] A further improvement or preferred embodiment of the aforementioned method for solving the gear pitting dynamics model based on microscopic surface topography, wherein the roughness is the arithmetic mean difference of the profile of the measured surface, which is defined as: ; where yi is the height of each roughness peak, measured by a profilometer, and n refers to the total number of roughness peaks.

[0100] A further improvement or preferred implementation of the aforementioned method for solving the gear pitting dynamics model based on microscopic surface morphology is that in step (4b), the solution method used is the fourth-order Runge-Kutta method. BRIEF DESCRIPTION OF THE DRAWINGS

[0101] Figure 1 This is a flow chart of a method for solving the gear pitting dynamics model based on microscopic surface morphology;

[0102] Figure 2 It is a schematic diagram of a gear meshing pair with fractal characteristics;

[0103] Figure 3 This is a physical picture of the optical profilometer;

[0104] Figure 4 It is the actual surface image of the gear taken by the optical profilometer during profile analysis;

[0105] Figure 5 It is an optical profilometer that extracts three-dimensional roughness maps;

[0106] Figure 6 is the optical profilometer to extract the 2D roughness map (0 hours);

[0107] Figure 7 is an optical profilometer to extract a two-dimensional roughness map (400 hours);

[0108] Figure 8 is an optical profilometer to extract a two-dimensional roughness map (600 hours);

[0109] Figure 9 is the result of calculation of single-side tooth backlash;

[0110] Figure 10 is the comparison of contact stiffness of gear pair after correction;

[0111] Figure 11 It is the cantilever beam model of the gear tooth when the base circle radius is smaller than the root circle radius;

[0112] Figure 12 It is a cantilever beam model when the base circle radius is larger than the root circle radius;

[0113] Figure 13 It is a cantilever beam model of internal gear;

[0114] Figure 14 It refers to the schematic diagram of lubricated contact of gear pairs considering surface roughness;

[0115] Figure 15 It is a representation of the thickness of the lubricating oil film. DETAILED DESCRIPTION

[0116] The present invention is described in detail below with reference to specific embodiments.

[0117] One of the prerequisites for dynamic characteristic analysis of a gear transmission system with rough tooth surfaces is to establish an accurate dynamic model. An accurate and reasonable dynamic model can simulate the gear working process in a comprehensive and detailed manner, thereby streamlining calculations while improving the accuracy of the dynamic response of the gear transmission system.

[0118] like Figure 2 As shown in the figure, the gear surface obtained by machining will have a processing texture, which is not smooth at the microscopic level and has many rough profile peaks. At the same time, as the gear working time increases, the gear roughness peak profile will become more complex. The height of these rough profile peaks will directly affect the tooth side clearance and normal contact stiffness of the meshing gears, and will also affect the friction, wear and lubrication conditions of the gears, and ultimately affect the dynamic characteristics of the gears. Therefore, in order to obtain more accurate results, some nonlinear internal excitation and hybrid elastohydrodynamic lubrication need to consider the influence of tooth surface roughness.

[0119] To this end, this application establishes a tooth surface morphology and pitting defect model based on the following steps, namely:

[0120] Step 1: Establish tooth surface morphology and pitting defect model

[0121] 1a. First, sample the selected gear tooth surface morphology and analyze the nonlinear excitation factors inside the gear transmission system. The sampling method of the gear tooth surface morphology is usually optical measurement. In this embodiment, the μ-scan optical profiler of Mangophy Optical Technology Co., Ltd. is used. Figure 3 As shown in the figure, it uses multi-aperture confocal technology, replacing the detector with a rotating disk with multiple apertures. During operation, the objective lens moves vertically, using a technique similar to tomography, to measure and acquire three-dimensional surface data of a workpiece in a very short time (approximately a few seconds). This optical profilometer is widely used to measure workpiece surface topography and conduct three-dimensional topography research at the micron or even nanometer scale, such as three-dimensional surface topography, two-dimensional depth topography, profile (depth, width, curvature, angle), and surface roughness. Parameter information is shown in Table 1.

[0122] Table 1 Optical profiler parameter information

[0123]

[0124] In this example, the sun gear used in the gear transmission system of a model wind turbine speed increaser was selected for solution verification. Gear surface topography can vary significantly under different operating conditions. In this example, the sun gear was measured after operating for 0 hours, 400 hours, and 600 hours. The gear to be analyzed has 27 teeth, a module of 12, and is precision ground.

[0125] Figure 4 The gear surface was photographed during profile analysis using an optical profilometer. It can be seen that the surface of the unoperated gear is relatively flat, and only clear grinding marks can be seen due to grinding. After 400 hours of operation, more meshing marks appear on the tooth surface, and slight wear occurs in some places. After 600 hours of operation, the tooth surface shows more severe wear than after 400 hours of operation, and more micro-pitting pits appear on the tooth surface.

[0126] The roughness components extracted from the three tooth surfaces are as follows Figure 5 As shown in the figure, it shows that when the gear enters the stable wear stage after the running-in stage, the surface of the part will be damaged, the degree of wear will increase, and the surface of the part will become increasingly rough.

[0127] The true two-dimensional roughness is extracted from the three-dimensional roughness. Figures 6-8 As shown in the figure, the Ra value is the arithmetic mean difference of the profile of the measured surface, which is defined as: ;

[0128] where y i is the height of each roughness peak, measured by the profilometer, Refers to the total number of rough peaks.

[0129] 1b. Establish a mathematical model of gear tooth surface morphology based on fractal dimension D and fractal roughness G

[0130] The fractal dimension D determines the frequency content of the surface contour waviness and the fractal roughness G represents a microscopic height scale parameter that is independent of frequency. The above two fractal parameters are related to the power spectrum. The relationship can be expressed as: ;

[0131] Where ω is the spatial frequency, γ is the constant related to the spatial frequency, and γ is generally taken as 1.5. The fractal dimension and fractal roughness under different roughness are shown in Table 2:

[0132] Table 2 Fractal dimension and fractal roughness at different roughness

[0133]

[0134] The calculation methods of the fractal dimension D and the fractal roughness G are as follows:

[0135] Taking the logarithm of the spatial frequency lgω as the abscissa and the logarithm of the power spectral density lgP(ω) as the ordinate to plot a graph, obtaining the corresponding slope and intercept, and thus calculating the corresponding fractal dimension D and fractal roughness G;

[0136] Fractal theory is a method for describing rough surfaces with fractal characteristics. It has a certain self-similar structure and can be used to characterize the rough surfaces in gear pairs. Different from the roughness obtained by statistical methods, the fractal method has the advantage of not depending on the instrument resolution. In fractal theory, the fractal dimension D and the fractal roughness G are the main parameters for determining the surface topography. The smoother the part surface, the larger the fractal dimension D and the smaller the fractal roughness G. According to fractal theory, the rough meshing interface has fractal characteristics and can be simulated by a series of micro-protrusions, such as Figure 2 the micro-protrusions in

[0137] At present, the methods for characterizing rough surfaces using fractal theory mainly include the fractal Brown function simulation method, the fractal interpolation simulation method, the W-M function simulation method, the inverse Fourier transform simulation method, the time series simulation method, and the composite fractal simulation method, etc. Among them, the W-M function simulation method is the most commonly used.

[0138] 1c. Establish a mathematical model of gear tooth surface deformation

[0139] As shown in Figure 2 in and respectively represent the fractal rough surfaces of the driving gear and the driven gear, represents the clearance between the two tooth profiles without considering the surface topography. Since most gear pairs need to be lubricated through lubricants, etc., and there will be certain errors in the processing and assembly of gears, there is always a backlash in the gear meshing pair. By measuring the surface topography of the gear, it can be found that it has fractal characteristics, so it can also be expressed by the W-M function. Then the expression of the fractal rough surface is: ;

[0140] In the formula, z i (t), i = p or g, are the fractal rough surfaces of the two gears; λ is the scale coefficient, D is the fractal dimension, where 1 < D < 2, and n is the frequency exponent, which determines the length dimension of the micro-protrusions; the above formula does not take the fractal roughness G into account. To obtain a more accurate result, the above formula is modified to:

[0141] ;

[0142] In the formula, L g is the sampling length of the fractal rough surface. For meshing gears, Lg is the actual contact length, expressed as ;

[0143] Among them, R m =R p R g / (R p +R g ) is the comprehensive curvature radius of the meshing pair, R p is the curvature radius of the driving wheel, R g F is the curvature radius of the driven wheel, N is the meshing force;

[0144] Then the single-side gear backlash considering the surface topography can be obtained as:

[0145]

[0146] Where, is the single-sided tooth backlash of the j-th gear pair, is the initial single-sided tooth clearance without considering the surface topography.

[0147] Take the tooth side clearance =25μm is 12.5 μm, then the calculated single-side tooth side clearance is as follows Figure 9 As shown in the figure, it can be seen that with the increase of working time, the tooth surface becomes rougher and rougher, the fractal dimension D value becomes smaller, the fractal roughness G value becomes larger, and the unilateral tooth clearance of the gear pair changes more. When Ra=0.96, the maximum fluctuation amplitude can reach 4μm.

[0148] When the absolute value of the meshing displacement is less than the tooth side clearance, the gear teeth will not deform, so the meshing displacement and meshing force are not in a linear relationship. Let δj represent the projection of the relative displacement of the j-th gear pair on the meshing line, then the deformation between the gear pairs can be expressed by the piecewise function express:

[0149]

[0150] The meshing force between gear pair j is It can be expressed as: ;in Represents the meshing stiffness of the jth gear pair.

[0151] Step 2: Establish a time-varying mesh stiffness correction model

[0152] 2a. Establish a mathematical model for optimizing contact stiffness considering surface topography

[0153] Mesh stiffness is the ability of meshing teeth to resist deformation due to meshing forces. It varies with the number of meshing tooth pairs. Since the contact ratio of involute spur gears is typically between 1 and 2, the meshing process is a process of alternating single and double tooth engagements. Therefore, the mesh stiffness of a gear pair is also time-varying. When the mesh stiffness of a gear is time-varying, it can lead to significant nonlinearity in the gear system. Therefore, this step is primarily used to model and correct this time-varying mesh stiffness.

[0154] The gear teeth are regarded as nonlinear variable cross-section cantilever beams, which will produce deformation during the meshing process. There are several basic potential energies in the meshing process, including the Hertz contact potential energy U h , bending potential energy U b , shear potential energy U s , radial compression potential energy U a , and the corresponding stiffnesses are Hertz contact stiffness k h , bending stiffness k b , shear stiffness k s , radial compression stiffness k a ;

[0155] When considering the rough surface, the Hertz contact stiffness is mainly affected by the surface micromorphology. Therefore, this application mainly focuses on the correction of the Hertz contact stiffness.

[0156] Based on the traditional algorithm, assuming that the gear teeth are isotropic elastic bodies, the Hertz contact stiffness expression can be obtained as: ; Where E is the equivalent elastic modulus, ν is the Poisson's ratio, and L is the tooth width;

[0157] It can be seen that in the traditional calculation method, the Hertzian contact stiffness is only related to the gear material and tooth width. However, since the tooth surface is rough, the line contact of the cylinder is also uneven at the microscopic level. The meshing interface is not a smooth contact but a mutual contact of rough surfaces. That is, the traditional algorithm ignores the influence of surface morphology on contact stiffness and cannot accurately describe the meshing characteristics of gears with different rough surfaces.

[0158] On the basis of step 1, based on fractal theory, at the microscopic scale, the stiffness of the contact interface is the combination of the stiffness of all asperities. During the elastic deformation stage of the asperities, the stiffness of a single asperity can be expressed as: ; Where α is the actual contact area of ​​a single micro-convex body;

[0159] For a single rough body, elastic deformation occurs first at the initial contact, followed by plastic deformation. The contact area divided into two stages is expressed as: Where σ y is the yield strength, l is the length of the asperity;

[0160] The critical deformation of the two-stage division can be expressed as: ;

[0161] The height of a single asperity can be expressed as: ;

[0162] When the critical deformation exceeds the height of the micro-convex body, the contact will be in the pure elastic deformation stage, and the critical length of the micro-convex body l c It can be expressed as: ;

[0163] On a pair of contacting tooth surfaces, the distribution function of the asperities is: ; In the formula, α represents the contact area of ​​a single micro-convex body, α L represents the maximum contact area of ​​a single asperity, λ m represents the surface coefficient;

[0164] The expression is as follows:

[0165] ;

[0166] ;

[0167] Where, L a is the diameter of the bottom surface of the microconvex body

[0168] Combining the above formulas, based on the contact area and critical length, the contact stiffness considering the surface topography can be obtained by using quadratic integration:

[0169] ;

[0170] Where D and G are the equivalent fractal dimension and equivalent fractal roughness of the meshing pair, respectively, and ; ;

[0171] Where D p With D g are the fractal dimensions of the driving wheel and the driven wheel in the gear pair, G p With G g are the fractal roughness of the driving wheel and the driven wheel in the gear pair respectively;

[0172] Substituting the data measured by the profiler and the gear data into the formula, the contact stiffness of each selected gear can be obtained. In this embodiment, the contact stiffness is as follows: Figure 10 shown.

[0173] It can be seen from this that the contact stiffness is greatly affected by the surface. For different fractal dimensions and fractal roughness, the contact stiffness changes significantly. As the tooth surface becomes rougher, the profile arithmetic mean difference Ra increases, the fractal dimension D decreases, and the fractal roughness G increases, the gear tooth contact stiffness decreases from 8.6708×109N / m to 4.69139×109N / m.

[0174] 2b. Establish a mathematical model for optimizing meshing stiffness considering surface topography

[0175] According to the variable cross-section cantilever beam and the involute tooth profile characteristics of the gear, the bending stiffness k can be derived b , shear stiffness k s , radial compression stiffness k a However, in actual implementation, the gear transmission system contains both external and internal meshing, and the calculation of bending stiffness, shear stiffness and radial compression stiffness of the two are also different, so they need to be analyzed separately.

[0176] (2B1) Mathematical model when the base circle diameter is smaller than the root circle radius

[0177] like Figure 11 The figure shows the cantilever beam model of the gear tooth when the base circle radius is smaller than the root circle radius. At this time, there is no transition curve in the gear tooth. a With F b The meshing force F Nj The components of force in the horizontal and vertical directions; d(y) represents the distance between the meshing point and the base circle in the direction of tooth height; h x represents the distance from the base circle x to the center line of the gear tooth; h(y) represents the distance between the meshing point and the center line of the gear tooth in the direction of tooth thickness; h f Indicates the distance from the intersection of the root circle and the tooth profile to the center line of the gear tooth; h b Indicates the distance from the starting point of the involute on the base circle to the center line of the gear tooth; Indicates the meshing force F Nj With F b Angle; Indicates the angle between the starting point of the involute on the base circle and the center line and the center line of the gear tooth; Indicates the angle between the starting point of the involute on the root circle and the center line of the gear tooth; Indicates the angle between the tooth centerline and the line connecting the base circle tangent points; Indicates the pitch circle pressure angle, represents the root circle pressure angle; r b Indicates the radius of the gear base circle; r f Indicates the tooth root circle radius;

[0178] Then the bending potential energy Ub can be expressed as: ;

[0179] The shear potential energy Us can be expressed as: ;

[0180] The radial compression potential energy Ua can be expressed as: ;

[0181] Where, E and G represent the gear elastic modulus and shear modulus respectively, I x A represents the moment of inertia of the tooth section at a distance x from the base circle in the tooth height direction. x represents the cross-sectional area of ​​the above section;

[0182] Combine Figure 11 The parameter derivation formula is as follows:

[0183] ;

[0184] ;

[0185] ;

[0186] ;

[0187] ;

[0188] ;

[0189] Substituting the parameter derivation formula into the potential energy expression, the relationship between each stiffness and parameter can be obtained as follows:

[0190] ;

[0191]

[0192] ;

[0193] ;

[0194] In the above formula, the positive and negative signs of the upper and lower limits of the integral are defined as follows: Figure 11 As shown, when the angle is below the center line, its value is defined as negative, and when it is above the center line, it is defined as positive.

[0195] When the base circle radius is smaller than the root circle radius, the meshing stiffness is calculated using the above method. The gear tooth deformation between the root circle and the base circle that does not exist will be taken into account, which will cause the calculated meshing stiffness to be smaller than the actual value. In order to obtain accurate meshing stiffness, the integral upper limit of the original formula needs to be changed from Change to , its expression is as follows,

[0196]

[0197]

[0198]

[0199] (2B2) Mathematical model when the base circle radius is greater than the root circle radius

[0200] like Figure 12 As shown, when the base circle radius is larger than the root circle radius, compared with Figure 11 , a transition curve appears between the root circle and the base circle. If the simplified algorithm is used, the length of the cantilever beam of the gear teeth in the model will be reduced, causing the calculated value of the gear tooth meshing stiffness to be larger than the actual value. Therefore, in order to improve the calculation accuracy, the gear tooth deformation between the root circle and the base circle needs to be taken into account in the calculation model. Then:

[0201]

[0202]

[0203]

[0204]

[0205]

[0206]

[0207] Refers to the distance between the pitch circle and the tooth root position. It refers to the distance from a point on the transition curve from tooth root to involute to the center line of the gear tooth. Indicates the tooth root transition curve radius;

[0208] Further calculation yields:

[0209] ;

[0210] in ;

[0211] ;

[0212] ;

[0213] (2B3) Mathematical model of internal gear

[0214] The tooth groove of the internal meshing gear is involute, so the cantilever beam model uses the tooth groove center line as the auxiliary calculation line, such as Figure 13 As shown in the figure, D ris the root of the inner gear ring, F is the meshing point, X is D r Any point on the line connecting points X and F, x represents the distance between points X and F in the direction of the center line, h represents the distance between points X and F in the direction of the center line, x is the distance from point X to the center line of the inner ring gear teeth, d is the distance between the meshing point and the tooth root in the center line direction, It is the angle between the tangent line of the base circle through point X and the center line of the tooth groove. The angle between the tangent line of the base circle at point F and the center line of the tooth groove is, for The angle between its component perpendicular to the center line, , according to the angle relationship in the cantilever beam, the corresponding calculation method is transformed to:

[0215] ;

[0216] ;

[0217] ;

[0218] in

[0219]

[0220]

[0221]

[0222]

[0223] It refers to the base circle radius; refers to the tooth root radius; B refers to the deformation coefficient; It refers to the angle between the starting point of the involute on the root circle and the center line of the gear tooth;

[0224] Step 3: Establish a mathematical model for the coupled analysis of hybrid elastohydrodynamic lubrication of pitting gears under rough surfaces

[0225] In actual operation, since lubricating oil is added between gears to reduce friction and wear, an oil film is formed between the gear teeth. The oil film has a certain oil film stiffness and oil film damping during operation, which makes the meshing stiffness and meshing damping different from the stiffness and damping derived in Chapter 2 without considering lubrication. At the same time, the friction coefficient will also be affected to a certain extent after lubrication, thereby affecting the gear friction. Therefore, it is very necessary to perform fluid analysis on the lubricating oil. Considering the surface roughness, the lubricated contact of the spur gear pair is as follows: Figure 14As shown in the figure, Z1(x) and Z2(x) represent the surface morphology of the meshing teeth respectively. It can be found that the lubrication condition of the lubricating oil flowing through the two surfaces is greatly affected by the surface roughness. Therefore, this step is mainly based on the relationship between the curvature radius, dynamic load, oil film entrainment speed, etc. under different surface roughness and the change of meshing position to establish a mathematical model for the hybrid elastohydrodynamic lubrication coupling analysis of pitting gears under rough surfaces. The mathematical model is based on the Reynolds equation and includes:

[0226] (3a) The lubrication model of the gear system under stable working condition can be expressed as:

[0227] ;

[0228] Where x represents the coordinate in the same direction as the main speed, ρ represents the density of the lubricating oil, h represents the thickness of the lubricating oil film, η represents the viscosity of the lubricating oil, p represents the oil film pressure on the lubricating oil, and u o1 、u o2 Represent the velocities of the upper and lower surfaces along the x direction, u o Indicates the entrainment speed;

[0229] The physical meaning of this model is the balance of flow. The left side of the model is the pressure flow (Poiseuille flow) formed by the pressure gradient of the lubricating oil film along the x direction, and the rightmost side of the model is the shear flow (Couette flow) caused by the surface velocity.

[0230] (3b) Mathematical model of film thickness in mixed elastohydrodynamic lubrication process

[0231] like Figure 15 As shown in the figure, during the mixed elastohydrodynamic lubrication process, the tooth surfaces of the meshing gears will squeeze each other, so the elastic deformation of the tooth surfaces needs to be considered in the calculation process. Therefore, the surface elastic deformation needs to be taken into account in the film thickness equation, and the mathematical model is established as follows:

[0232]

[0233] Where h(x) is the lubricating oil film thickness at x, h0 is the center film thickness, R is the instantaneous curvature radius, υ(x) is the elastic deformation at x, and s r (x) is the surface roughness at x. The occurrence of pitting directly determines the change in roughness. Therefore, the change in gear meshing lubrication characteristics can be indirectly affected by the change in micro-roughness.

[0234] (3b) Mathematical model of viscosity and density in mixed elastohydrodynamic lubrication process

[0235] The viscosity of lubricating oil is one of the key physical properties. It can represent the flow resistance encountered by gears during operation. During operation, the viscosity of lubricating oil is not only related to the lubricating oil itself, but also to the temperature and pressure to which it is subjected. When the working pressure of common mineral oil exceeds 20MPa, the viscosity begins to change significantly with pressure, and the rate of change of viscosity also increases with increasing pressure. The viscosity-pressure-temperature equation is derived from experimental results, and the mathematical model can be expressed as:

[0236]

[0237] Where T0 is the initial temperature, is the lubricant viscosity at T=T0, is the pressure viscosity coefficient, generally taken as =1.96×108, is a constant, usually taken =0.68;

[0238] The density of lubricating oil is also greatly affected by pressure. It increases with increasing pressure. At the same time, the increase in temperature will also cause the thermal expansion of the lubricating oil, resulting in an increase in volume and a decrease in density. The mathematical model of density changes with pressure and temperature can be expressed as:

[0239]

[0240] Where, is the lubricant density when p=0, T=T0, is the density-temperature coefficient, generally taken ;

[0241] Step 4: Solve the nonlinear dynamic mathematical model of gears with rough surfaces

[0242] Surface roughness has a significant impact on the time-varying meshing stiffness, tooth side clearance, and mixed elastohydrodynamic lubrication of gears. At the same time, the oil film generated when lubrication is introduced into the gear transmission system will also affect dynamic internal excitations such as friction, meshing stiffness, and meshing damping. The combined effect of these dynamic internal excitations will make the dynamic characteristics of the time-varying system exhibit very strong nonlinearity, and also make the movement of the gear transmission system more complex. This step establishes a solution model for the dynamic characteristics of the gear transmission system with rough tooth surfaces. Based on the solution model, combined with the Poincare section and power spectrum, the nonlinear dynamic mathematical model of the gear with rough surfaces is solved;

[0243] Based on the combination of steps 1 to 3, the dynamic characteristics solution process of the gear transmission considering the pitting of the rough tooth surface is as follows: Figure 1 As shown, the specific steps include:

[0244] (4a) The gear tooth roughness was measured using an optical profilometer, and the basic parameters of the gear and the gear excitation frequency were prepared. The tooth side clearance considering the surface topography, the time-varying mesh stiffness considering the surface topography, the comprehensive error excitation, and the tooth surface friction torque were calculated.

[0245] (4b) Substitute the calculated dynamic internal excitation into the derived dynamic differential equations and solve them using the fourth-order Runge-Kutta method to obtain the relative meshing displacement of the system. At the beginning of the iteration, the friction torque is set to 0.

[0246] (4c) After calculating the relative meshing displacement and relative meshing speed, the relative entrainment speed, meshing force, and equivalent curvature radius are determined based on the gear parameters, relative meshing displacement, and relative meshing speed. The kinematic and load parameters are substituted into the hybrid elastohydrodynamic lubrication model every Δt to iterate and obtain the hybrid elastohydrodynamic lubrication characteristics at that moment. The variation pattern of the hybrid elastohydrodynamic lubrication characteristics within a meshing cycle is calculated.

[0247] (4d) Substitute the obtained time-varying hybrid elastohydrodynamic lubrication characteristics into the dynamic mathematical model for the next iteration. The iterative convergence condition is:

[0248] ;

[0249] Where Q is the number of iterations, and γ is taken as all meshing positions calculated within one meshing cycle;

[0250] is the gear meshing force, It refers to time, Refers to the gear that bears the force. Refers to the number of gears.

[0251] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit the scope of protection of the present invention. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention may be modified or replaced by equivalents without departing from the essence and scope of the technical solutions of the present invention.

Claims

1. A method for solving the gear pitting dynamics model based on microscopic surface morphology, characterized in that: It includes steps 1 to 4. The specific contents of each step include: Step 1. Establish the tooth surface morphology and pitting defect model; specifically divided into the following sub-steps: 1a. Sample the selected gear tooth surface morphology and collect the nonlinear excitation parameters inside the gear transmission system; 1b. Establish a mathematical model of gear tooth surface morphology based on fractal dimension D and fractal roughness G, which can be expressed as: ; in, is the power spectrum, ω is the spatial frequency, and γ is a constant with the spatial frequency; The calculation method of fractal dimension D and fractal roughness G is as follows: The logarithm of spatial frequency lgω is used as the horizontal coordinate, and the logarithm of power spectrum density lgP(ω) is used as the vertical coordinate to obtain the corresponding slope and intercept, and the corresponding fractal dimension D and fractal roughness G are calculated; 1c. Establish a mathematical model of gear tooth surface deformation The WM function is combined with fractal roughness to express the gear tooth surface deformation, and the expression of the fractal rough surface is obtained as follows: ; where z i (t), i = p or g, respectively refer to the fractal rough surfaces of the driving wheel and the driven wheel; λ is the scale coefficient, D is the fractal dimension, 1 < D < 2, n is the frequency exponent, refers to the maximum frequency exponent; and respectively represent the fractal rough surfaces of the driving wheel and the driven wheel, represents the clearance between two tooth profiles without considering the surface topography, refers to the scale coefficient, and t refers to time; L g is the fractal rough surface sampling length, expressed as ; Among them, R m =R p R g / (R p +R g ) is the comprehensive curvature radius of the meshing pair, R p is the curvature radius of the driving wheel, R g F is the curvature radius of the driven wheel, N is the meshing force; E is the elastic modulus; The gear unilateral backlash considering the surface topography is further obtained as: Where, is the single-sided tooth backlash of the j-th gear pair, is the initial single-sided tooth backlash without considering the surface topography; It refers to the fractal dimension of the driving wheel and the driven wheel. It refers to the fractal dimension of the driving wheel and the driven wheel; use Represents the projection length of the relative displacement of the jth pair of gears on the meshing line, then the deformation between the gear pairs is expressed by the piecewise function Expressed as: Meshing force between gear pairs j Expressed as: ;in represents the meshing stiffness of the j-th gear pair; is the derivative of the projected length of the relative displacement of the jth gear pair on the meshing line; Step 2: Establish a time-varying mesh stiffness correction model 2a. Establish a mathematical solution model for contact stiffness optimization considering surface topography Consider the gear teeth as a nonlinear cantilever beam with variable cross-section. Based on step 1 and fractal theory, the stiffness of the contact interface is considered as the combination of the stiffness of all asperities. Based on the contact area and critical length, the contact stiffness considering the surface topography can be obtained by using quadratic integration. : ； Where D and G are the fractal dimension and fractal roughness of meshing, respectively, and ; ; represents the stiffness of a single asperity and ; is the actual contact area of ​​a single asperity; is the critical length of the asperity and , is Poisson's ratio, is the yield strength; is the distribution function of the asperity, Represents the maximum contact area of ​​a single micro-convex body, La is the diameter of the bottom surface of the micro-convex body, Represents the surface coefficient; D p is the fractal dimension of the driving wheel in the gear pair, D g is the fractal dimension of the driven wheel, G p With G g are the fractal roughness of the driving wheel and the driven wheel in the gear pair respectively; 2b. Establish a mathematical model for optimizing meshing stiffness considering surface topography, including: (2B1) Mathematical model when the base circle diameter is smaller than the root circle radius In the mathematical model when the base circle diameter is smaller than the root circle radius, define F a With F b The meshing force F Nj The components of force in the horizontal and vertical directions; d(y) represents the distance between the meshing point and the base circle in the direction of tooth height; h x represents the distance from the base circle x to the center line of the gear tooth; h(y) represents the distance between the meshing point and the center line of the gear tooth in the direction of tooth thickness; h f Indicates the distance from the intersection of the root circle and the tooth profile to the center line of the gear tooth; h b Indicates the distance from the starting point of the involute on the base circle to the center line of the gear tooth; Indicates the meshing force F Nj With F b The angle between Indicates the angle between the starting point of the involute on the base circle and the center line and the gear tooth center line; Indicates the angle between the starting point of the involute on the root circle and the center line of the gear tooth; Indicates the angle between the tooth centerline and the line connecting the base circle tangent points; Indicates the pitch circle pressure angle, represents the root circle pressure angle; r b Indicates the radius of the gear base circle; r f Indicates the tooth root circle radius; Then the bending potential energy Ub can be expressed as: ; The shear potential energy Us can be expressed as: ; The radial compression potential energy Ua can be expressed as: ; Where, E, Represent the gear elastic modulus and shear modulus, I x A represents the moment of inertia of the tooth section at a distance x from the base circle in the tooth height direction. x represents the cross-sectional area of ​​the above section; The parameter derivation formula is as follows: ； ； ； ； ； ； Substituting the parameter derivation formula into the potential energy expression, we can obtain the Hertz contact stiffness , bending stiffness , shear stiffness The relationship with parameters: ； ； ； ； ； ； (2B2) Mathematical model when the base circle radius is greater than the root circle radius When the base circle radius is larger than the root circle radius, the gear tooth deformation between the root circle and the base circle is taken into account in the calculation model, and the result is: ； ； ； ； Refers to the distance between the pitch circle and the tooth root position. It refers to the distance from a point on the transition curve from tooth root to involute to the center line of the gear tooth. Indicates the tooth root transition curve radius; Further we get: ； in ; ； ； (2B3) Mathematical model of internal gear In the mathematical model of internal gears, define D r is the root of the inner gear ring, F is the meshing point, X is D r Any point on the line connecting points X and F, x represents the distance between points X and F in the direction of the center line, h represents the distance between points X and F in the direction of the center line, x is the distance from point X to the center line of the inner ring gear teeth, d is the distance between the meshing point and the tooth root in the center line direction, It is the angle between the tangent line of the base circle through point X and the center line of the tooth groove. The angle between the tangent line of the base circle at point F and the center line of the tooth groove is, for The angle between the component perpendicular to the center line and the component perpendicular to the center line is calculated by transforming the corresponding calculation method according to the angle relationship in the cantilever beam: ； ； ； in , It refers to the base circle radius; refers to the tooth root radius; B refers to the deformation coefficient; It refers to the angle between the starting point of the involute on the root circle and the center line of the gear tooth; Step 3: Establish a mathematical model for the coupled analysis of hybrid elastohydrodynamic lubrication of pitting gears under rough surfaces Based on the relationship between the curvature radius, dynamic load, oil film entrainment speed and other factors as the meshing position changes under different surface roughness, a mathematical model for the hybrid elastohydrodynamic lubrication coupling analysis of pitting gears under rough surfaces is established, including: (3a) The lubrication model of the gear system under stable working condition is expressed as: ； In the lubrication model of the gear system under stable working conditions, x represents the coordinate in the same direction as the main speed, ρ represents the density of the lubricating oil, h represents the thickness of the lubricating oil film, η represents the viscosity of the lubricating oil, p represents the oil film pressure on the lubricating oil, and u o1 、u o2 Represent the velocities of the upper and lower surfaces along the x direction, u o represents the entrainment velocity. The physical meaning of the lubrication model under the stable working state of the gear system is the balance of flow. The left side of the model is the pressure flow formed by the pressure gradient of the lubricating oil film along the x direction, and the rightmost side of the model is the shear flow caused by the surface velocity. (3b) The mathematical model of film thickness in the mixed elastohydrodynamic lubrication process is expressed as: In the film thickness mathematical model of the mixed elastohydrodynamic lubrication process, h(x) is the lubricating oil film thickness at x, h0 is the central film thickness, R is the instantaneous curvature radius, υ(x) is the elastic deformation at x, and s r (x) is the surface roughness at position x; (3b) Mathematical model of viscosity and density in mixed elastohydrodynamic lubrication process The viscosity mathematical model of the mixed elastohydrodynamic lubrication process is obtained by fitting the experimental results and is expressed as: In the viscosity mathematical model of the mixed elastohydrodynamic lubrication process, T0 is the initial temperature, is the lubricant viscosity at T=T0, is the pressure viscosity coefficient, is a constant; The density mathematical model of the mixed elastohydrodynamic lubrication process is expressed as: In the density mathematical model of the mixed elastohydrodynamic lubrication process, is the lubricant density when p=0, T=T0, is the density-temperature coefficient; Step 4: Solve the nonlinear dynamic mathematical model of gears with rough surfaces Based on steps 1 to 3, the solution process for the dynamic characteristics of the gear transmission with pitting corrosion considering the rough tooth surface includes: (4a) Use an optical profilometer to detect gear tooth roughness and confirm the basic parameters and nonlinear excitation parameters of the gear; (4b) Substituting the calculated dynamic internal excitation into the derived corresponding mathematical model for iterative solution to solve the relative meshing displacement and relative meshing speed of the system. At the beginning of the iteration, the friction torque is set to 0; (4c) After the relative meshing displacement and relative meshing speed are calculated, the relative entrainment speed, meshing force, and equivalent curvature radius are determined based on the basic parameters and nonlinear excitation parameters of the gear, the relative meshing displacement, and the relative meshing speed. The kinematic and load parameters are substituted into the hybrid elastohydrodynamic lubrication model every Δt to iterate and obtain the hybrid elastohydrodynamic lubrication characteristics at that moment. The hybrid elastohydrodynamic lubrication characteristic data within one meshing cycle are calculated. (4d) Substitute the obtained time-varying hybrid elastohydrodynamic lubrication characteristic data into the dynamic mathematical model for the next iteration until the model converges. The iterative convergence condition is: ； Where Q is the number of iterations, and γ is taken as all meshing positions calculated within one meshing cycle; is the gear meshing force, It refers to time, Refers to the gear that bears the force. Refers to the number of gears.

2. The method for solving the gear pitting dynamics model based on microscopic surface morphology according to claim 1 is characterized in that: The nonlinear excitation parameters in step (1a) include three-dimensional surface topography, two-dimensional depth topography, profile (depth, width, curvature, angle), and surface roughness.

3. The method for solving the gear pitting dynamics model based on microscopic surface morphology according to claim 2 is characterized in that: The roughness is the arithmetic mean difference of the profile of the measured surface, which is defined as: ; where yi is the height of each roughness peak, measured by the profilometer, Refers to the total number of rough peaks.

4. The method for solving the gear pitting dynamics model based on microscopic surface morphology according to claim 1 is characterized in that: In the step (4b), the solution method used is the fourth-order Runge-Kutta method.

Images