Quasi-periodic state-dependent regenerative milling stability analysis method and system
By constructing state-dependent time-delay differential equations and multi-scale dynamic analysis, the analytical challenge of quasi-periodic vibration in flexible systems in existing technologies has been solved, enabling accurate prediction of milling stability and optimization of process parameters.
Patent Information
- Application Number
- CN202610527720.6
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2026-04-21
- Publication Date
- 2026-08-25
AI Technical Summary
Existing milling stability analysis methods cannot accurately handle quasi-periodic vibrations in flexible systems, resulting in stability lobe diagrams that cannot fully cover actual machining scenarios and cannot provide reliable process parameter guidance for high-speed CNC milling.
A state-dependent time-delay differential equation model is constructed, and a multi-scale dynamic analysis framework is adopted. By solving the model in layers at the principal spindle period and large period scale, the characteristic parameters of quasi-periodic vibration are identified, the milling stability limit diagram is drawn, and the basis for optimizing process parameters is provided.
It achieves unified modeling and stability determination of periodic chatter and quasi-periodic chatter, improves the prediction accuracy of chatter critical depth of cut and boundary shape, ensures computational efficiency and self-consistency of results, and is applicable to high-speed CNC milling of flexible systems.
Smart Images

Figure CN122634827A_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of CNC milling technology, and more specifically, relates to a quasi-periodic state-dependent regenerative milling stability analysis method and system. Background Technology
[0002] With the widespread adoption of high-speed CNC milling technology in high-end manufacturing fields such as aerospace, large wind tunnel components, and molds, chatter is becoming increasingly prevalent in tool-workpiece systems operating at high speeds and large cutting depths. Chatter is a self-excited vibration caused by the coupling of multiple factors, including structural flexibility and the regenerative effect of cutting forces. It not only leads to deterioration of part surface roughness and deviations from design accuracy, but also significantly shortens tool life and, in severe cases, can even cause tool breakage or machine tool instability, becoming a key bottleneck restricting the improvement of machining quality and production efficiency. To mitigate chatter risks, engineering practice urgently needs to define the stable cutting range under different machining parameters. The Stability Lobe Diagram (SLD), as a core tool, can visually present the maximum stable depth of cut corresponding to different spindle speeds, providing important guidance for optimizing process parameters. Traditional stability analysis is mostly based on the periodic motion assumption that "vibration is synchronized with spindle speed". This assumption has certain applicability in systems with high rigidity, but in systems with greater flexibility, the actual vibration often presents a quasi-periodic pattern of "high-frequency flutter superimposed with low-frequency envelope". That is, the self-excited flutter frequency and the spindle speed are not simply integer multiples of each other. This poses a serious challenge to traditional analysis methods and also creates an urgent need for more accurate stability analysis technology.
[0003] Currently, the core theoretical foundation for milling chatter stability analysis is regenerative cutting theory. This theory posits that the instantaneous cutting thickness is jointly influenced by the current feed motion and the "ripple" surface left by the previous machining. Consequently, a time delay term related to the tooth passage cycle is introduced into the dynamic equation. Various technical solutions have been developed around this theory, resulting in related patents. For example, patent document CN119830569A discloses a milling chatter stability prediction method based on an improved fully discrete method. This method establishes a dynamic differential equation containing regenerative effects and transforms it into a state-space equation. Third-order Hermite interpolation and second-order Newton interpolation are used to approximate the state term and time delay term, respectively. Then, a Floquet transition matrix is constructed to determine system stability, thus improving computational accuracy and efficiency. Patent document CN118502355A, targeting micro-milling scenarios, proposes a stability simulation prediction method considering process damping and tool runout. A feedback dynamic model is established based on the undeformed cutting thickness, and stability prediction is achieved through a fully discrete method using a direct integration scheme. From the perspective of technical path classification, existing technologies mainly include three categories: First, the frequency domain / analytical method, which transforms the regeneration term into frequency domain parameters through equivalent transfer function derivation, and solves the limit cutting depth by combining the Nyquist criterion, etc. It has a fast calculation speed and is suitable for quickly drawing SLD diagrams; Second, the time domain-discretization method, which discretizes the spindle period into sub-steps and constructs a single-period discrete operator, and judges the stability of the periodic solution by using Floquet multipliers. It can take into account the direction coefficient matrix of periodic changes and multimodal coupling, and is widely used in engineering software; Third, the improved model that considers the periodic state-dependent effect, which introduces state-dependent time delay to fit the actual cutting trajectory, but still does not escape the constraints of the periodic assumption.
[0004] Although existing technologies have achieved certain application effects in specific scenarios, they still have significant shortcomings: (1) Frequency domain / analytical methods rely on the assumption of constant time delay, which cannot adapt to the actual working conditions where the time delay changes dynamically with the vibration state; (2) Both the time domain-discrete method and the periodic state-dependent improved model take "vibration and spindle speed synchronization" as the core premise, and can only handle periodic flutter where the self-excited frequency and the spindle frequency are integer multiples of each other. They lack effective analysis methods for quasi-periodic vibrations that are common in flexible systems, and it is difficult to accurately capture the stability boundary under such working conditions; (3) Although some improved models attempt to introduce state-dependent characteristics, due to the immaturity of analytical solutions and stability analysis frameworks, they can only be limited to the processing of periodic excitation flutter and cannot realize the system analysis of quasi-periodic flutter.
[0005] These shortcomings limit the stability lobe diagrams generated by existing methods, making it difficult to fully cover complex vibration scenarios in actual machining and providing reliable process parameter guidance for high-speed milling of flexible systems. This invention addresses these deficiencies by proposing an analysis method based on state-dependent time delay and large-period single-frequency harmonic balance. This method eliminates the need for pre-assuming low-frequency modulation frequencies and automatically identifies characteristic parameters of quasi-periodic vibrations, achieving unified modeling and stability determination of periodic and quasi-periodic flutter, thus filling a gap in existing technologies for quasi-periodic flutter analysis. Summary of the Invention
[0006] To address the aforementioned deficiencies or improvement needs of existing technologies, this invention provides a quasi-periodic state-dependent regenerative milling stability analysis method and system. Basic parameters are input as the foundational data for modeling. First, a state-dependent time-delay differential equation model suitable for flexible milling systems is constructed. Then, a multi-scale dynamic analysis framework is established, encompassing the basic vibration trajectory at the spindle periodic scale and the quasi-periodic vibration trajectory at the large periodic scale. The quasi-periodic trajectory further includes self-excited vibration components that are not integer multiples of the spindle frequency. Next, the multi-scale framework is solved hierarchically: first, using one spindle revolution as a benchmark, the periodic equilibrium solution is represented by finite harmonic series, and the periodic stable vibration trajectory of the spindle is solved through numerical iteration and residual convergence. Then, the Floquet theory is used to analyze the... The eigenvalues of the trajectory are used to identify the self-excited frequencies that are not integer multiples of the spindle frequency. The large period of the quasi-periodic motion is determined through rational approximation. The explicit harmonic balance method is used to discretize the quasi-periodic vibration at the large period scale, constructing an extended vector containing the current state and multiple delayed states to simulate the entire vibration response process under the coupling of the self-excited frequency and the spindle frequency. Finally, the stability analysis results of the quasi-periodic vibration are output. The results include at least the spindle periodic balance solution, the quasi-periodic balance solution, and the spectral radius of the large-period single-period discrete operator. Based on the results, the stability state of the milling system under different spindle speed-depth of cut parameters is determined, and a milling stability limit diagram containing the quasi-periodic chatter boundary is drawn, providing a basis for optimizing high-speed CNC milling process parameters.
[0007] To achieve the above objectives, according to a first aspect of the present invention, a quasi-periodic state-dependent method for analyzing the stability of regenerative milling is provided, comprising the following steps: S100: Establish state-dependent time-delay differential equations applicable to milling systems, input basic parameters as basic data for modeling, and establish unit models for simulation verification; S200: Based on one spindle revolution cycle, the periodic equilibrium solution is represented by finite harmonic series. The harmonic coefficients are solved by numerical iteration, and the self-excited frequency that is not an integer multiple of the spindle frequency is identified by combining Floquet theory. S300: Calculate the ratio of the self-excited frequency to the principal axis angular frequency, make a rational approximation within the preset maximum denominator range, and determine the large period of the quasi-periodic motion; S400: Using the large period as the fundamental frequency, the displacement response is expressed as a finite harmonic series. The state-dependent time delay and cutting thickness at each time node within the large period are solved. The fixed point condition is established through iteration to obtain the quasi-periodic equilibrium solution. S500: Linearize the system near the large-cycle equilibrium solution, construct an extended vector containing the current state and multi-level delayed states, construct single-step linear discrete operators step by step and multiply them to obtain a large-cycle single-cycle discrete operator; S600: Based on the spectral radius of the large-period single-period discrete operator, the stability of each working point is determined, and the milling stability limit diagram (SLD) containing the quasi-periodic unstable boundary is drawn on the spindle speed-depth of cut parameter plane.
[0008] Further, step S100 includes: S101: Input the basic data for modeling, including cutting force coefficient, climb and conventional milling settings, number of teeth, entry angle, exit angle, tool radius and other model parameters; S102: The tool-workpiece system is simplified into a mechanical system with vibrational degrees of freedom in both the tangential and normal directions, and a two-degree-of-freedom vibration equation is established: S103: Establish the geometric relationship between the tool path and the workpiece surface, and define the spatial constraints between the current position of each tooth and the historical regeneration trace; S104: Derive the state-dependent time-delay constraint relationship; S105: Determine the instantaneous cutting thickness based on the regeneration geometry: S106: Establish the cutting force model and calculate the tool orientation coefficient; S107: Output state-dependent time-delay differential equation: Integrating the vibration equation, cutting force, and time-delay constraint, the final dynamic equation is obtained as follows: ,in For vibration displacement, These are the mass, damping, and stiffness matrices, respectively, and t is time.
[0009] Furthermore, in step S101, the dynamic equation for the absence of cutting is: ; The state-dependent time-delay constraint relationship is as follows: ,in Indicates the vibration state. Represents an unknown, state-dependent time delay that changes with the state. The expression for instantaneous cutting thickness is: ; The linear cutting force model is ,in, k 0 Indicates the cutting edge force coefficient. k 1 The direction coefficient matrix represents the cutting force coefficient, which is determined by the instantaneous entry and exit angles of the cutting teeth and the tool geometry.
[0010] Further, step S200 includes: S201: Input the spindle speed and depth of cut, and introduce the state-dependent milling dynamics equation established in step S100; S202: Set the spindle cycle and select the harmonic cutoff order: Import the spindle speed, cutting depth, and the previously established state-dependent vibration equation; S203: Construct a harmonic expansion with the principal axis period frequency as the fundamental frequency, set the principal axis period and harmonic truncation order to control the calculation accuracy, and use the principal axis angular frequency as the fundamental frequency to express the vibration response as a harmonic series. S204: Discrete time nodes over the entire spindle cycle; S205: Solve for the state-dependent time delay at each time node based on geometric constraints; S206: Calculate the instantaneous cutting thickness of each tooth and select a physically reasonable time delay-thickness cutting combination; S207: Calculate the cutting force within the cycle; S208: Substitute the harmonic expansion of the cutting force and the structure into the dynamic equation to form the harmonic balance residual, and determine whether the residual converges. S209: If the residual converges, proceed to the self-excited frequency identification algorithm module; if the residual does not converge, return to S203. The self-excited frequency identification algorithm process includes: S210: Based on the equilibrium solution, linearize the state-dependent time-delay dynamic equation and construct the principal axis periodic discrete operator; S211: Calculate the eigenvalues of the discrete operator, i.e., the Floquet multiplier spectrum; S212: Identify complex conjugate pairs with membrane values close to 1 from the Floquet multiplier spectrum and extract the self-excited frequency; S213: Outputs the spindle periodic equilibrium solution and estimates the self-excited frequency.
[0011] Furthermore, the principal axis period is used as the reference period, and the periodic equilibrium solution is represented by a finite harmonic series: ,in a 0 For static offset, k is the harmonic order. Main axis angular frequency, H To cut off the harmonic number, t is time. Let be the Fourier coefficients to be determined.
[0012] Furthermore, by substituting the harmonic expansion into the aforementioned dynamic equation and treating the time delay and cutting thickness as functions of the harmonic coefficients, a closed-loop mapping of "harmonic coefficients → trajectory → time delay and cutting thickness → cutting force → trajectory derivative" is obtained. By comparing the closed-loop mapping with the trajectory derivative obtained analytically from harmonic expansion, the fixed-point conditions for the periodic equilibrium solution can be constructed. , where is the vector of all harmonic coefficients. As a residual operator, the harmonic coefficients corresponding to the principal axis periodic equilibrium solution can be obtained through numerical fixed-point iteration. .
[0013] Further, in step S300, the calculation formula for the large period is as follows: ,in Let q be the period of one spindle revolution, and q be the integer denominator obtained by rational approximation, satisfying the ratio of the self-excited frequency to the spindle angular frequency. .
[0014] Furthermore, step S400 utilizes the large period determined in S300 to transform the original quasi-period into a period based on... For the fundamental frequency of the large-cycle periodic motion, a new equilibrium solution is obtained within the large cycle. The fundamental frequency of the large cycle is defined as: ; The finite harmonic series expression for the displacement response is: ,in For the fundamental frequency, To cut off the harmonic number, a 0 This is the static offset, where t is time. Let be the harmonic coefficient to be determined; The fixed-point condition under the large periodic scale is established as follows: ,in Given the large-period harmonic coefficient vector, an equilibrium solution over the large period can be obtained through iterative solving. .
[0015] Furthermore, in step S500, the system variables are expressed near the large-period equilibrium solution as follows: A first-order Taylor expansion is performed on the cutting force term and the state-dependent geometric constraints to obtain a linear time-varying system with respect to displacement perturbation and historical perturbation. The single-step linear discrete operator includes the evolution of the current master state, the shift update of each level of delayed state, and the effect of delay stiffness interpolation.
[0016] According to a second aspect of the present invention, a quasi-periodic state-dependent regenerative milling stability analysis system is provided, comprising: The dynamics modeling module is used to establish state-dependent time-delay differential equations, and input basic parameters as the basic data for modeling and simulation verification. The self-excited frequency identification module is used to solve the spindle periodic equilibrium solution and identifies the self-excited frequency that is not an integer multiple of the spindle frequency based on Floquet theory. The large period determination module is used to make a rational approximation of the ratio of the self-excited frequency to the main shaft frequency to determine the large period of the quasi-periodic motion. The quasi-periodic solution module is used to solve quasi-periodic equilibrium solutions on a large-periodic scale, establish large-periodic fixed-point conditions, and iteratively converge. The discrete operator construction module is used to linearize the system near the large-period equilibrium solution, construct large-period single-period discrete operators and compute Floquet multipliers; The stability determination and SLD drawing module is used to determine the stability of each working point based on the operator spectrum radius and draw the milling stability limit diagram containing the quasi-periodic unstable boundary.
[0017] In summary, compared with the prior art, the above-described technical solutions conceived by this invention can achieve the following beneficial effects: 1. The method of the present invention does not require prior assumption of low-frequency modulation frequency. It automatically extracts the self-excited flutter frequency through harmonic balance and Floquet characteristic analysis of the principal axis single-cycle scale. Combined with rational approximation relationship, it constructs large-cycle quasi-periodic solution, which can accurately determine the characteristic period and frequency of self-excited quasi-periodic vibration. 2. The computational framework of this invention can simultaneously handle periodic chatter and quasi-periodic chatter of the main shaft. When the self-excitation frequency is an integer ratio to the main shaft frequency, it automatically degenerates into traditional periodic chatter analysis. When it is a non-integer rational ratio, it gives the large-period quasi-periodic solution and its stability boundary, so that the stability lobe diagram can fully reflect the two types of chatter and has a stronger ability to interpret low-frequency vibration phenomena. 3. The method of the present invention establishes a state-dependent time-delay constraint equation, explicitly considers the real geometric relationship between the cutting teeth and the historical regenerated surface, selects a physically reasonable cutting thickness and time-delay combination by adopting the minimum cutting thickness principle, and obtains the equivalent direction coefficient and stiffness increment that change with time after linearization. Compared with the traditional fixed time-delay or simplified geometric model, it can more accurately reflect the influence of factors such as entry / exit, tool runout and multi-tooth meshing, and improve the prediction accuracy of chatter critical depth of cut and boundary shape. 4. The method of this invention adopts frequency domain solution and linearization analysis uniformly on both single-cycle and large-cycle scales of the main shaft, and iteratively corrects the relationship between the self-excited frequency and the main shaft frequency. Under the premise of ensuring computational efficiency, it achieves the self-consistency and numerical robustness of the stability prediction results, and has good engineering applicability and promotion value. Attached Figure Description
[0018] Figure 1 This is a flowchart illustrating the overall algorithm workflow of the quasi-periodic state-dependent regenerative milling stability analysis method in this embodiment of the invention. Figure 2 This is a flowchart of state-dependent milling dynamics modeling in an embodiment of the present invention; Figure 3 This is a flowchart of the principal axis periodic equilibrium solution and Floquet self-excited frequency analysis in an embodiment of the present invention. Detailed Implementation
[0019] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention. Furthermore, the technical features involved in the various embodiments of this invention described below can be combined with each other as long as they do not conflict with each other.
[0020] like Figure 1 As shown, this embodiment takes the high-speed milling of a certain aerospace component as an example, and uses the quasi-periodic state-dependent regenerative milling stability analysis method of the present invention to perform stability analysis. The specific steps are as follows: S100: Establish a state-dependent milling dynamics model; establish state-dependent time-delay differential equations that can reflect regeneration effects, multi-tooth meshing of the tool, and trajectory-dependent geometric constraints, providing a unified dynamic basis for subsequent equilibrium solutions and linearization analysis.
[0021] S200: Analyze the periodic equilibrium solution of the main shaft and the Floquet self-excited frequency; solve the periodic equilibrium solution over one rotation of the main shaft, and identify the self-excited frequency related to regenerative flutter using Floquet theory, providing initial frequency information for large-period construction.
[0022] S300: Rational approximation of self-excited frequency and determination of large period; The ratio of self-excited frequency to principal axis frequency obtained in S200 is rationally approximated to determine the basic large period of quasi-periodic motion on a longer time scale, thereby providing a time structure for establishing a "single-frequency periodic model on a large period".
[0023] Self-excitation frequency estimation obtained from the previous module Calculate the ratio between its frequency and the principal axis angular frequency: In the pre-set maximum denominator Within the scope, with For the target value, the search satisfies: Integer pairs And prioritize rational number approximations with the smallest error and denominators not exceeding the maximum range. Within this framework, This can be viewed as a frequency-locked relationship between quasi-periodic vibration and spindle rotation. From the physical meaning, we can see that when the system simultaneously has the principal axis frequency and self-excitation frequency At that time, its overall motion in time It exhibits periodicity, that is, in After the principal axis of the circle rotates, the phase relationship of the system returns to its original state. Therefore, it can be... Defined as the large cycle of quasi-periodic motion.
[0024] The rational approximation parameters output by this module and the corresponding major cycles This is the core basis for large-period harmonic balance modeling and quasi-periodic stability analysis.
[0025] S400: Solve for the large-cycle single-frequency harmonic equilibrium quasi-periodic solution; using the large cycle determined in S300, transform the original quasi-periodic solution into a solution with... The fundamental frequency is a large-cycle periodic motion, in the large cycle Solve for a new equilibrium solution internally.
[0026] Define the fundamental frequency of a large period as: The displacement response is expressed as... Finite harmonic series of the fundamental frequency: in Let be the harmonic coefficients to be determined. Due to the increased time base, some components that previously exhibited "quasi-periodic" characteristics become integer harmonics at the new fundamental frequency, thus allowing for uniform processing in the frequency domain. Within the time interval... At the discrete time nodes, combined with the same state-dependent geometric constraints, the corresponding state-dependent time delay and cutting thickness are solved for each tooth and each time node, resulting in a time delay-cutting thickness pair covering the entire large cycle. The difference from the spindle cycle stage is that the historical tool position needs to span [a specific time point]. The phase relationship of historical trajectories in geometric constraints is more complex, but it can still be calculated consistently in the frequency domain through a unified large-period harmonic expansion.
[0027] Substituting the large-cycle harmonic expansion into the state-dependent time-delay differential equation, the cutting force within the large cycle is calculated using the obtained time delay and cutting thickness sequence. Substituting this back into the dynamic equation yields the displacement and velocity derivatives within the large cycle. Simultaneously, the time derivatives of displacement and velocity can be analytically obtained using the harmonic expansion. Comparing these two derivatives allows for the establishment of new fixed-point conditions on the large-cycle scale. in This represents the large-period harmonic coefficient vector. Through iterative solving, an equilibrium solution over the large period can be obtained. In the original physical time sense, it corresponds to a stable orbit of a quasi-periodic motion.
[0028] S500: Large-period linearization and construction of single-period discrete operators: Linearization is performed near the large-period equilibrium solution, and a single-period discrete operator is constructed to describe the perturbation evolution within the entire large period, thereby realizing the Floquet stability analysis of the quasi-periodic equilibrium solution.
[0029] Large-cycle equilibrium solution Nearby, write the system variable as A first-order Taylor expansion of the cutting force term and state-dependent geometric constraints yields the following results: and historical disturbances Linear time-varying systems, where The time delay is the time delay under equilibrium solution. At this point, the direction coefficient matrix and equivalent stiffness periodically change with time. Changes. Due to state-dependent time delays, historical disturbances need to be discretized in time into several levels of delayed states with fixed step sizes to form a finite-dimensional extended state vector. By constraining and reducing the time delay range, the number of delay levels can be guaranteed to be finite, thus making the monodromy matrix dimension controllable. A uniform time step is used throughout the entire large cycle. Within each time step, a linear discrete operator is constructed based on the linearized stiffness matrix and equivalent direction coefficients of the current and previous moments. This operator includes: 1. The evolution of the current master state in the next moment; 2. The "shifting" and updating of each level of delayed state; 3. The interpolation effect from the delay stiffness (reflected through weighting coefficients).
[0030] Multiplying the single-step operators at each time step in chronological order yields a single-period discrete operator describing the evolution of a disturbance from its initial state to the end state within a large period. The eigenvalues of this matrix are the Floquet multipliers of the quasi-periodic equilibrium solution. The large-period monodromy matrix and the corresponding Floquet multiplier spectrum output by this module are direct evidence for judging the stability of the quasi-periodic equilibrium solution.
[0031] S600: Determine quasi-periodic stability and draw stability leaf diagram: On the plane with given spindle speed and cutting depth parameters, repeat the above module process to determine the stability type at each working point, and finally form a stability leaf diagram that takes into account both periodic and quasi-periodic instability boundaries.
[0032] Specifically, a two-dimensional parametric mesh is constructed within a predetermined spindle speed range and depth of cut range. For each parameter point in the mesh, the following steps are performed in sequence: 1. State-dependent dynamics modeling and time-scale initialization; 2. Solving the spindle periodic equilibrium solution and identifying the Floquet self-excited frequency; 3. Rational approximation of the self-excited frequency and determination of the large period; 4. Solving the large-period single-frequency harmonic equilibrium quasi-periodic equilibrium solution; 5. Linearization of the large period and construction of the monodromy matrix; 6. Floquet multiplier calculation and stability determination.
[0033] Based on the spectral radius of the large-period monodromy matrix :when When the quasi-periodic equilibrium solution is considered stable under this parameter, it corresponds to the stable cutting zone; when At this time, the system is considered to be chattering unstable. Mapping the stable / unstable state and chatter type of each parameter point to the spindle speed–depth of cut plane yields a stability lobe diagram containing quasi-periodic instability boundaries.
[0034] like Figure 2 As shown, S100 is the core component of the connection state-dependent dynamic modeling and quasi-periodic analysis, enabling accurate solution of the principal axis periodic equilibrium and extraction of key quasi-periodic frequencies. Specific steps include: S101: Input basic parameters: Provides basic modeling data including structural dynamic parameters, such as mass matrix M, damping matrix C, stiffness matrix K, etc., as well as tool and machining geometry parameters, such as number of teeth, tool radius, cutting parameters, etc. S102: Establishing a two-degree-of-freedom vibration equation: The tool-workpiece system is simplified into a mechanical system with vibrational degrees of freedom in both the tangential and normal directions. Its dynamics without cutting can be written as: in For vibration displacement, These are the mass, damping, and stiffness matrices, respectively.
[0035] S103: Establish the geometric relationship between the tool path and the workpiece surface, and define the spatial constraints between the current position of each tooth and the historical regeneration trace; S104: State-dependent time-delay constraint relationship: During milling, the actual cutting thickness of each tooth is determined by the current tool position, the historical tool positions of the previous few revolutions, and the feed motion. Since the location of the regenerated mark is related to the tool vibration state, the time interval experienced by the cutting teeth passing through the same point is no longer a constant, but an unknown, state-dependent time delay that varies with the state. For each cutting tooth, the geometric constraints can be written as follows: This constraint simultaneously characterizes the relationship between the tool path, workpiece contour, and spindle rotation angle, thus implicitly defining... With vibration state The function is continuous.
[0036] S105: Determine the instantaneous cutting thickness based on the regeneration geometry. , usually S106: Establish the cutting force model and calculate the tool direction coefficient: Based on S105, a linear cutting force model is adopted, for example... The direction coefficient matrix is determined by the instantaneous entry and exit angles of the cutting teeth and the tool geometry. This cutting force is then substituted back into the dynamic equations to form a state-dependent time-delay differential equation (SDDDE).
[0037] S107: Output state-dependent time-delay differential equation: Integrating the vibration equation, cutting force, and time-delay constraint, the final dynamic equation is obtained as follows: like Figure 3 The diagram shows the flowchart for "Spindle Period Equilibrium Solution and Self-Excited Frequency Identification" in milling chatter analysis. The core task is to solve for the stable vibration trajectory within the spindle period and identify the key frequencies of quasi-periodic chatter. It includes two modules: "Spindle Period Equilibrium Solution Algorithm" and "Self-Excited Frequency Identification Algorithm." The "Spindle Period Equilibrium Solution Algorithm" flowchart includes: S201: Input spindle speed and depth of cut, and introduce the state-dependent milling dynamics equation established in the previous module; Select one spindle revolution cycle As a reference period, the periodic equilibrium solution Represented using finite harmonic series: in Main axis angular frequency, To cut off the harmonic number, Let be the Fourier coefficients to be determined.
[0038] S202: Set the spindle cycle and select the harmonic cutoff order: Import the spindle speed, cutting depth, and the previously established state-dependent vibration equation; S203: Construct a harmonic expansion with the principal axis period frequency as the fundamental frequency: Set the principal axis period and harmonic truncation order to control the calculation accuracy, and use the principal axis angular frequency as the fundamental frequency to express the vibration response as a harmonic series. S204: Discrete time nodes over the entire spindle cycle; S205: Solve for the state-dependent time delay at each time node based on geometric constraints; S206: Calculate the instantaneous cutting thickness of each tooth and select a physically reasonable time delay-thickness cutting combination; For each time node, the geometric constraint equations Considered as a time-delay variable The nonlinear algebraic equations are solved numerically or rewritten into first-order differential form and integrated numerically to obtain the reasonable time delay sequence of each tooth at each time node. Simultaneously, based on the obtained time delay and displacement response, the instantaneous cutting thickness of each tooth is calculated, and all candidate cutting thicknesses are screened to select the regenerated tooth sequence corresponding to the minimum positive cutting thickness, thereby obtaining a physically reasonable time delay-cutting thickness combination.
[0039] S207: Calculate the cutting force within the cycle; S208: Substitute the harmonic expansion of the cutting force and the structure into the dynamic equation to form the harmonic balance residual, and determine whether the residual converges. Substituting the harmonic expansion into the dynamic equations, and treating time delay and cutting thickness as functions of the harmonic coefficients, we obtain a closed-loop mapping: "harmonic coefficients → trajectory → time delay and cutting thickness → cutting force → trajectory derivative". Comparing this mapping with the trajectory derivative obtained analytically from the harmonic expansion, we can construct the fixed-point conditions for the periodic equilibrium solution. in Let be the vector of all harmonic coefficients. This is the residual operator. The harmonic coefficients corresponding to the principal axis periodic equilibrium solution can be obtained through numerical fixed-point iteration (such as iterative updates in the least squares sense). .
[0040] S209: If the residual converges, proceed to the self-excited frequency identification algorithm module; if the residual does not converge, return to S203. The self-excited frequency identification algorithm process includes: S210: Based on the equilibrium solution, linearize the state-dependent time-delay dynamic equation and construct the principal axis periodic discrete operator; In the obtained equilibrium solution Nearby, the cutting force and geometric constraints are linearized to the first order to obtain information about the displacement perturbation. For a linear time-varying system, the regeneration term is represented as a delay stiffness matrix with time-dependent engagement angles and trajectories. To handle the time delay term, the historical states are discretized and expanded into a finite-dimensional "delayed state sequence" at fixed time steps, constructing an extended state vector containing the current state and several levels of delayed states. Within one principal axis period, by progressively constructing and multiplying single-step linear operators, a single-period discrete operator (monodromy matrix) describing the perturbation evolution within one period can be obtained. .
[0041] S211: Calculate the eigenvalues of the discrete operator, i.e., the Floquet multiplier spectrum; S212: Identify complex conjugate pairs with membrane values close to 1 from the Floquet multiplier spectrum and extract the self-excited frequency; S213: Outputs the spindle periodic equilibrium solution and estimates the self-excited frequency.
[0042] right Find the eigenvalues to obtain the Floquet multiplier set. If a pair of conjugate complex multipliers Its magnitude is close to 1, and its phase angle corresponding frequency This is the candidate value of the self-excitation frequency added to the equilibrium solution of the principal axis period. If the multiplier magnitude is significantly greater than 1, it indicates that the system has undergone regenerative flutter on the principal axis period scale. If the magnitude is close to 1 but deviates from the principal axis and its harmonic frequencies, it can be determined that there is a self-excitation frequency that is "not an integer multiple of the principal axis period", which provides a basis for subsequent quasi-periodic analysis.
[0043] Those skilled in the art will readily understand that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.
Claims
1. A quasi-periodic state-dependent method for analyzing the stability of regenerative milling, characterized in that, Includes the following steps: S100: Establish state-dependent time-delay differential equations applicable to milling systems, input basic parameters as basic data for modeling, and establish unit models for simulation verification; S200: Based on one spindle revolution cycle, the periodic equilibrium solution is represented by finite harmonic series. The harmonic coefficients are solved by numerical iteration, and the self-excited frequency that is not an integer multiple of the spindle frequency is identified by combining Floquet theory. S300: Calculate the ratio of the self-excited frequency to the principal axis angular frequency, make a rational approximation within the preset maximum denominator range, and determine the large period of the quasi-periodic motion; S400: Using the large period as the fundamental frequency, the displacement response is expressed as a finite harmonic series. The state-dependent time delay and cutting thickness at each time node within the large period are solved. The fixed point condition is established through iteration to obtain the quasi-periodic equilibrium solution. S500: Linearize the system near the large-cycle equilibrium solution, construct an extended vector containing the current state and multi-level delayed states, construct single-step linear discrete operators step by step and multiply them to obtain a large-cycle single-cycle discrete operator; S600: Based on the spectral radius of the large-period single-period discrete operator, the stability of each working point is determined, and the milling stability limit diagram (SLD) containing the quasi-periodic unstable boundary is drawn on the spindle speed-depth of cut parameter plane.
2. The quasi-periodic state-dependent regenerative milling stability analysis method according to claim 1, characterized in that, Step S100 includes: S101: Input the basic data for modeling, including cutting force coefficient, climb and conventional milling settings, number of teeth, entry angle, exit angle, tool radius and other model parameters; S102: The tool-workpiece system is simplified into a mechanical system with vibrational degrees of freedom in both the tangential and normal directions, and a two-degree-of-freedom vibration equation is established: S103: Establish the geometric relationship between the tool path and the workpiece surface, and define the spatial constraints between the current position of each tooth and the historical regeneration trace; S104: Derive the state-dependent time-delay constraint relationship; S105: Determine the instantaneous cutting thickness based on the regeneration geometry: S106: Establish the cutting force model and calculate the tool orientation coefficient; S107: Output state-dependent time-delay differential equation: Integrating the vibration equation, cutting force, and time-delay constraint, the final dynamic equation is obtained as follows: ,in For vibration displacement, These are the mass, damping, and stiffness matrices, respectively, and t is time.
3. The quasi-periodic state-dependent regenerative milling stability analysis method according to claim 2, characterized in that, In step S101, the dynamic equation for the absence of cutting is: ; The state-dependent time-delay constraint relationship is as follows: ,in Indicates the vibration state. Represents an unknown, state-dependent time delay that changes with the state. The expression for instantaneous cutting thickness is: ; The linear cutting force model is ,in, k 0 Indicates the cutting edge force coefficient. k 1 The direction coefficient matrix represents the cutting force coefficient, which is determined by the instantaneous entry and exit angles of the cutting teeth and the tool geometry.
4. The quasi-periodic state-dependent regenerative milling stability analysis method according to any one of claims 1-3, characterized in that, Step S200 includes: S201: Input the spindle speed and depth of cut, and introduce the state-dependent time-delay differential equation; S202: Set the spindle cycle and select the harmonic cutoff order: Import the spindle speed, cutting depth, and the previously established state-dependent time-delay differential equation; S203: Construct a harmonic expansion with the principal axis period frequency as the fundamental frequency, set the principal axis period and harmonic truncation order to control the calculation accuracy, and use the principal axis angular frequency as the fundamental frequency to express the vibration response as a harmonic series. S204: Discrete time nodes over the entire spindle cycle; S205: Solve for the state-dependent time delay at each time node based on geometric constraints; S206: Calculate the instantaneous cutting thickness of each tooth and select a physically reasonable time delay-thickness cutting combination; S207: Calculate the cutting force within the cycle; S208: Substitute the harmonic expansion of the cutting force and the structure into the dynamic equation to form the harmonic balance residual, and determine whether the residual converges. S209: If the residual converges, proceed to the self-excited frequency identification algorithm module; if the residual does not converge, return to S203. The self-excited frequency identification algorithm process includes: S210: Based on the equilibrium solution, linearize the state-dependent time-delay dynamic equation and construct the principal axis periodic discrete operator; S211: Calculate the eigenvalues of the discrete operator, i.e., the Floquet multiplier spectrum; S212: Identify complex conjugate pairs with membrane values close to 1 from the Floquet multiplier spectrum and extract the self-excited frequency; S213: Outputs the spindle periodic equilibrium solution and estimates the self-excited frequency.
5. The quasi-periodic state-dependent regenerative milling stability analysis method according to claim 4, characterized in that, The principal axis period is used as the reference period, and the periodic equilibrium solution is represented by a finite harmonic series: ,in a 0 For static offset, k is the harmonic order. Main axis angular frequency, H To cut off the harmonic number, t is time. Let be the Fourier coefficients to be determined.
6. The quasi-periodic state-dependent regenerative milling stability analysis method according to claim 5, characterized in that, Substituting the harmonic expansion into the aforementioned dynamic equation, and considering the time delay and cutting thickness as functions of the harmonic coefficients, we obtain a closed-loop mapping of "harmonic coefficients → trajectory → time delay and cutting thickness → cutting force → trajectory derivative". By comparing the closed-loop mapping with the trajectory derivative obtained analytically from harmonic expansion, the fixed-point conditions for the periodic equilibrium solution can be constructed. , where is the vector of all harmonic coefficients. As a residual operator, the harmonic coefficients corresponding to the principal axis periodic equilibrium solution can be obtained through numerical fixed-point iteration. .
7. The quasi-periodic state-dependent regenerative milling stability analysis method according to any one of claims 1-3, characterized in that, In step S300, the calculation formula for the large period is as follows: ,in Let q be the period of one spindle revolution, and q be the integer denominator obtained by rational approximation, satisfying the ratio of the self-excited frequency to the spindle angular frequency. .
8. The quasi-periodic state-dependent regenerative milling stability analysis method according to claim 7, characterized in that, Step S400 utilizes the large period determined in S300 to transform the original quasi-period into a period based on the large period determined in S300. For the fundamental frequency of the large-cycle periodic motion, a new equilibrium solution is obtained within the large cycle. The fundamental frequency of the large cycle is defined as: ; The finite harmonic series expression for the displacement response is: ,in For the fundamental frequency, To cut off the harmonic number, a 0 This is the static offset, where t is time. Let be the harmonic coefficient to be determined; The fixed-point condition under the large periodic scale is established as follows: ,in Given the large-period harmonic coefficient vector, an equilibrium solution over the large period can be obtained through iterative solving. .
9. The quasi-periodic state-dependent regenerative milling stability analysis method according to claim 8, characterized in that, In step S500, the system variables are expressed as near the large-period equilibrium solution. A first-order Taylor expansion is performed on the cutting force term and the state-dependent geometric constraints to obtain a linear time-varying system with respect to displacement perturbation and historical perturbation. The single-step linear discrete operator includes the evolution of the current master state, the shift update of each level of delayed state, and the effect of delay stiffness interpolation.
10. A quasi-periodic state-dependent regenerative milling stability analysis system, characterized in that, include: The dynamics modeling module is used to establish state-dependent time-delay differential equations, and input basic parameters as the basic data for modeling and simulation verification. The self-excited frequency identification module is used to solve the spindle periodic equilibrium solution and identifies the self-excited frequency that is not an integer multiple of the spindle frequency based on Floquet theory. The large period determination module is used to make a rational approximation of the ratio of the self-excited frequency to the main shaft frequency to determine the large period of the quasi-periodic motion. The quasi-periodic solution module is used to solve quasi-periodic equilibrium solutions on a large-periodic scale, establish large-periodic fixed-point conditions, and iteratively converge. The discrete operator construction module is used to linearize the system near the large-period equilibrium solution, construct large-period single-period discrete operators and compute Floquet multipliers; The stability determination and SLD drawing module is used to determine the stability of each working point based on the operator spectrum radius and draw the milling stability limit diagram containing the quasi-periodic unstable boundary.
Citation Information
Patent Citations
Micro-milling machining stability simulation prediction method considering flutter influence in machining process
CN118502355A
Milling chatter stability prediction method
CN119830569A
