CN103823787B

CN103823787B - Multi-turning-tool parallel turning stability judgment method based on differential quadrature method

Info

Publication number: CN103823787B
Application number: CN201410060539.6A
Authority: CN
Inventors: 丁烨; 牛金波; 朱利民; 丁汉
Original assignee: Shanghai Jiao Tong University
Current assignee: Wuxi Liman Robot Technology Co ltd
Priority date: 2014-02-21
Filing date: 2014-02-21
Publication date: 2017-01-18
Anticipated expiration: 2034-02-21
Also published as: CN103823787A

Abstract

The multi-turn tool parallel turning stability determination method based on differential quadrature method provided by the present invention includes the steps of: performing dynamic modeling on the multi-turn tool parallel turning processing system, establishing a multi-time-delay second-order differential equation; establishing and obtaining normalization The state-space equations are optimized; on the adjacent unit intervals [0,1] and [‑1,0], the second kind of Chebyshev points are used as discrete points; using the differential quadrature method, based on the Lagrangian interpolation function , use the displacement item at the discrete point to represent the velocity item; determine the interval where the discrete point of the time delay item is located, and use the second type of Chebyshev point in the interval to represent the time delay item; construct the distance between the two adjacent unit intervals State transition matrix, according to the Floquet theory to determine the stability of the original system. Compared with the traditional single-tool turning process, the present invention adopts the differential quadrature method to analyze the dynamic characteristics of the multi-tool parallel turning system, obtains optimized cutting parameters, and greatly improves the machining efficiency.

Description

Stability Judgment Method for Multi-tool Parallel Turning Based on Differential Quadrature Method

技术领域technical field

本发明涉及多车刀并行车削稳定性判定的一种新方法，即微分求积法(Differential Quadrature Method)，具体为利用微分求积法合理选择切削参数使多车刀并行车削避免再生颤振的影响进行高效率高质量加工。The present invention relates to a new method for judging the stability of multi-tool parallel turning, that is, differential quadrature method (Differential Quadrature Method). Affects high-efficiency and high-quality processing.

背景技术Background technique

在机械加工领域，车削是最为常见和常用的加工方式之一。传统的单车刀加工研究已经非常成熟，而多把车刀并行车削的加工方式是近几年学术界才提出的新概念，尤其是多车刀并行车削稳定性的研究，尚处于起步阶段。理论上，多车刀并行车削的加工效率要远远高于传统车削，但是由于多车刀并行车削的加工机理要比传统车削复杂得多。影响多车刀并行车削加工质量的主要因素是加工稳定性，因此研究多车刀并行车削的动力学机理，合理选择加工参数，可以有效地避免再生颤振的发生，从而实现加工过程的平稳运行，在保证加工质量的同时实现工件的高效车削。In the field of machining, turning is one of the most common and commonly used processing methods. The research on traditional single turning tool machining has been very mature, but the processing method of parallel turning with multiple turning tools is a new concept proposed by the academic circle in recent years, especially the research on the stability of parallel turning with multiple turning tools is still in its infancy. In theory, the processing efficiency of parallel turning with multiple turning tools is much higher than that of traditional turning, but the processing mechanism of parallel turning with multiple turning tools is much more complicated than that of traditional turning. The main factor affecting the quality of multi-turning tool parallel turning is machining stability. Therefore, studying the dynamic mechanism of multi-turning tool parallel turning and selecting the processing parameters reasonably can effectively avoid the occurrence of regenerative chatter and realize the smooth operation of the machining process. , to achieve efficient turning of the workpiece while ensuring the processing quality.

文献1“E.Budak,E.Ozturk,Dynamics and stability of parallelturningoperations.CIRP Annals—Manufacturing Technology 60(2011)383-386.”利用频率法对并行车削进行了稳定性分析，并加以实验验证。文献中的多车刀固定在不同的刀架上，刀具之间不存在耦合效应。其步骤如下：Document 1 "E. Budak, E. Ozturk, Dynamics and stability of parallel turning operations. CIRP Annals—Manufacturing Technology 60 (2011) 383-386." The frequency method is used to analyze the stability of parallel turning operations and verify it experimentally. The multi-turning tools in the literature are fixed on different tool holders, and there is no coupling effect between the tools. The steps are as follows:

(1)建立多车刀并行车削的动力学模型；(1) Establish a dynamic model of parallel turning with multiple turning tools;

(2)获得稳定边界处动力学方程；(2) Obtain the dynamic equation at the stable boundary;

(3)利用搜索法得到多车刀并行车削在切削参数空间的稳定区域；(3) Using the search method to obtain the stable region of multi-tool parallel turning in the cutting parameter space;

(4)进行时域仿真验证及切削实验验证。(4) Carry out time-domain simulation verification and cutting experiment verification.

文献2“E.Ozturk,E.Budak,Modeling dynamics of parallel turningoperations.Proceedings of 4^th CIRP International Conference on HighPerformance Cutting,2010.”讨论了两种多车刀并行车削的情况：多把车刀安装于不同刀架之上；多把车刀安装于同一刀架之上。这两种情况的分析计算方法类似，与文献1的步骤相同。Document 2 "E.Ozturk, E.Budak, Modeling dynamics of parallel turning operations. Proceedings of 4th CIRP International Conference on ^{HighPerformance} Cutting, 2010." Discussed two cases of parallel turning with multiple turning tools: multiple turning tools installed in different On the tool holder; multiple turning tools are installed on the same tool holder. The analysis and calculation methods of these two cases are similar, and the steps are the same as those in Document 1.

频域法步骤简洁，计算速度快，但是不利于考虑各种复杂工况条件。时域仿真法结果直观，但是计算速度慢，只能得到单一切削参数组合条件下的加工稳定性，不利于绘制加工参数空间的稳定性图谱。与上述方法相比，Bellman等在二十世纪七十年代提出的微分求积法具有计算速度快、结果精度高等优点，已经被广泛用于各个工程技术领域。多车刀并行车削的动力学方程是多时滞微分方程，将微分求积法做适当推广并用于多车刀并行车削的加工稳定性判定具有重要的生产价值和意义。The frequency domain method has simple steps and fast calculation speed, but it is not conducive to considering various complex working conditions. The results of the time-domain simulation method are intuitive, but the calculation speed is slow, and it can only obtain the machining stability under the condition of a single cutting parameter combination, which is not conducive to drawing the stability map of the machining parameter space. Compared with the above methods, the differential quadrature method proposed by Bellman et al. in the 1970s has the advantages of fast calculation speed and high result accuracy, and has been widely used in various engineering technology fields. The dynamic equation of multi-tool parallel turning is a multi-delay differential equation. It is of great production value and significance to properly promote the differential quadrature method and apply it to the processing stability judgment of multi-tool parallel turning.

发明内容Contents of the invention

针对现有技术中的缺陷，本发明的目的是依据多车刀并行车削的动力学特性，提供一种切削过程稳定性判定的新方法，即微分求积法，为多车刀并行车削加工参数的选择提供有效依据，在保证无再生颤振高质量加工的前提下获得尽量高的加工效率，产生良好的经济效益。Aiming at the deficiencies in the prior art, the purpose of the present invention is to provide a new method for judging the stability of the cutting process based on the dynamic characteristics of multi-turning tool parallel turning, that is, the differential quadrature method, which is a multi-turning tool parallel turning processing parameter The selection provides an effective basis to obtain the highest possible processing efficiency under the premise of ensuring high-quality processing without regenerative chatter, resulting in good economic benefits.

本发明提供了一种基于微分求积法的多车刀并行车削稳定性判定方法，包括如下步骤：The invention provides a method for judging the stability of multi-turning tool parallel turning based on differential quadrature method, comprising the following steps:

对多车刀并行车削加工系统进行动力学建模，建立多时滞二阶微分方程。进行状态空间变化，建立状态空间方程，之后进行归一化处理，得到归一化后的标准方程；The dynamic modeling of multi-turning tool parallel turning system is carried out, and the multi-time-delay second-order differential equation is established. Perform state space changes, establish state space equations, and then perform normalization processing to obtain normalized standard equations;

在相邻的单位区间[0,1]和[-1,0]上以第二类切比雪夫点(Chebyshev-Gauss-Lobatto Points)为离散点，将标准方程等价地离散为一组代数方程；On the adjacent unit intervals [0,1] and [-1,0], using the second kind of Chebyshev points (Chebyshev-Gauss-Lobatto Points) as discrete points, the standard equation is equivalently discretized into a set of algebraic equation;

利用微分求积法，基于拉格朗日插值函数，用离散点处的位移项表示速度项；判定时滞项离散点所处区间，基于拉格朗日插值函数，用所在区间的第二类切比雪夫点表示时滞项；Using the differential quadrature method, based on the Lagrangian interpolation function, the displacement item at the discrete point is used to represent the velocity item; to determine the interval of the discrete point of the delay item, based on the Lagrangian interpolation function, use the second type of the interval Chebyshev points represent time-delay terms;

构造相邻两个单位区间[0,1]和[-1,0]之间的状态转移矩阵，根据Floquet理论判定原系统的稳定性；若状态转移矩阵的所有特征值的模均小于1，则系统是稳定的；若状态转移矩阵的任一特征值的模大于1，则系统是不稳定的；由此可以得到多车刀并行车削系统在切削参数空间的稳定性谱图。Construct the state transition matrix between two adjacent unit intervals [0,1] and [-1,0], and judge the stability of the original system according to the Floquet theory; if the modulus of all eigenvalues of the state transition matrix are less than 1, Then the system is stable; if the modulus of any eigenvalue of the state transition matrix is greater than 1, the system is unstable; thus the stability spectrum of the multi-tool parallel turning system in the cutting parameter space can be obtained.

具体地，根据本发明提供的基于微分求积法的多车刀并行车削稳定性判定方法，包括如下步骤：Specifically, according to the differential quadrature method provided by the present invention, the method for judging the stability of parallel turning with multiple turning tools includes the following steps:

步骤1：建立多车刀并行车削动力学方程，对多车刀并行车削动力学方程进行整理，得多时滞二阶微分方程；Step 1: Establish multi-tool parallel turning dynamic equations, sort out multi-turn tool parallel turning dynamic equations, and multi-time-delay second-order differential equations;

步骤2：对所述多时滞二阶微分方程进行状态空间变换，得到状态空间方程；Step 2: performing state-space transformation on the multi-delay second-order differential equation to obtain a state-space equation;

步骤3：对所述状态空间方程进行归一化处理，得到标准形式的状态空间方程；Step 3: performing normalization processing on the state-space equation to obtain the state-space equation in standard form;

步骤4：对标准形式的状态空间方程进行周期离散，将其等价转化为一组代数方程作为状态空间方程表达式；Step 4: Periodically discretize the state-space equation in standard form, and convert it into a set of algebraic equations as the expression of the state-space equation;

步骤5：基于拉格朗日插值函数，以单位区间[0,1]上的第二类切比雪夫点为离散点，利用微分求积法，用位移项表示状态空间方程表达式中的导数项；Step 5: Based on the Lagrangian interpolation function, take the second kind of Chebyshev point on the unit interval [0,1] as a discrete point, and use the differential quadrature method to represent the derivative in the expression of the state space equation with the displacement term item;

步骤6：基于拉格朗日插值函数，对状态空间方程表达式中的时滞项进行所在区间判断，若时滞项的离散点属于区间[0,1]，则以[0,1]上的第二类切比雪夫点为离散点表示时滞项；若时滞项的离散点属于区间[-1,0]，则以[-1,0]上的第二类切比雪夫点为离散点表示时滞项；Step 6: Based on the Lagrangian interpolation function, judge the interval of the delay item in the expression of the state space equation. If the discrete point of the delay item belongs to the interval [0,1], then use [0,1] The second kind of Chebyshev points on [-1,0] are discrete points representing time-delay items; Discrete points represent delay items;

步骤7：构造相邻单位区间[0,1]和[-1,0]之间的状态转移矩阵，根据Floquet理论判定原系统的稳定性。Step 7: Construct the state transition matrix between adjacent unit intervals [0,1] and [-1,0], and judge the stability of the original system according to the Floquet theory.

优选地，所述步骤1中，所述多车刀并行车削动力学方程，如公式(1)所示：Preferably, in the step 1, the multi-turning tool parallel turning dynamics equation, as shown in formula (1):

$\begin{matrix} [\begin{matrix} {\overset{\cdot\cdot \cdot\cdot}{z z}}_{11} ((t t)) \\ {\overset{\cdot\cdot \cdot\cdot}{z z}}_{22} ((t t)) \end{matrix}] + + [\begin{matrix} 22 {ζ ζ}_{11} {ω ω}_{n no 11} & 00 \\ 00 & 22 {ζ ζ}_{22} {ω ω}_{n no 22} \end{matrix}] [\begin{matrix} {\overset{\cdot \cdot}{z z}}_{11} ((t t)) \\ {\overset{\cdot \cdot}{z z}}_{22} ((t t)) \end{matrix}] + + [\begin{matrix} {ω ω}_{n no 11}^{22} & 00 \\ 00 & {ω ω}_{n no 22}^{22} \end{matrix}] [\begin{matrix} {z z}_{11} ((t t)) \\ {z z}_{22} ((t t)) \end{matrix}] = = \\ {K K}_{f f} [\begin{matrix} \frac{{ω ω}_{n no 11}^{22}}{{k k}_{11}} {a a}_{11} ((- - {z z}_{11} ((t t)) + + {z z}_{22} ((t t - - \frac{τ τ}{22})))) \\ \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} {a a}_{11} ((- - {z z}_{22} ((t t)) + + {z z}_{11} ((t t - - \frac{τ τ}{22})))) + + \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} (({a a}_{22} - - {a a}_{11})) ((- - {z z}_{22} ((t t)) + + {z z}_{22} ((t t - - τ τ)))) \end{matrix}] \end{matrix} - - - - - - ((11))$

其中，下标1表示第一车刀，下标2表示第二车刀，z(t)表示动力学响应位移，ω_n表示固有圆频率，K_f表示进给方向的切削力系数，τ表示工件旋转周期，ζ表示相对阻尼比，k表示刚度系数，a₁表示第一车刀的切削深度，a₂表示第二车刀的切削深度，表示动力学加速度响应，表示动力学速度响应，t表示时间；Among them, subscript 1 indicates the first turning tool, subscript 2 indicates the second turning tool, z(t) indicates the dynamic response displacement, ω _n indicates the natural circular frequency, K _f indicates the cutting force coefficient in the feed direction, and τ indicates The workpiece rotation period, ζ represents the relative damping ratio, k represents the stiffness coefficient, a ₁ represents the cutting depth of the first turning tool, a ₂ represents the cutting depth of the second turning tool, represents the dynamic acceleration response, Indicates dynamic velocity response, t indicates time;

所述多时滞二阶微分方程，如公式(2)所示：The multi-time-delay second-order differential equation, as shown in formula (2):

$\begin{matrix} [\begin{matrix} {\overset{\cdot\cdot \cdot\cdot}{z z}}_{11} ((t t)) \\ {\overset{\cdot\cdot \cdot\cdot}{z z}}_{22} ((t t)) \end{matrix}] + + [\begin{matrix} 22 {ζ ζ}_{11} {ω ω}_{11 n no} & 00 \\ 00 & 22 {ζ ζ}_{22} {ω ω}_{n no 22} \end{matrix}] [\begin{matrix} {\overset{\cdot &Center Dot;}{z z}}_{11} ((t t)) \\ {\overset{\cdot &Center Dot;}{z z}}_{22} ((t t)) \end{matrix}] + + [\begin{matrix} {ω ω}_{n no 11}^{22} + + {K K}_{f f} \frac{{ω ω}_{n no 11}^{22}}{{k k}_{11}} {a a}_{11} & 00 \\ 00 & {ω ω}_{n no 22}^{22} + + {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} {a a}_{22} \end{matrix}] [\begin{matrix} {z z}_{11} ((t t)) \\ {z z}_{22} ((t t)) \end{matrix}] = = \\ [\begin{matrix} 00 & {K K}_{f f} \frac{{ω ω}_{n no 11}^{22}}{{k k}_{11}} {a a}_{11} \\ {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} {a a}_{11} & 00 \end{matrix}] [\begin{matrix} {z z}_{11} ((t t - - \frac{τ τ}{22})) \\ {z z}_{22} ((t t - - \frac{τ τ}{22})) \end{matrix}] + + [\begin{matrix} 00 & 00 \\ 00 & {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} (({a a}_{22} - - {a a}_{11})) \end{matrix}] [\begin{matrix} {z z}_{11} ((t t - - τ τ)) \\ {z z}_{22} ((t t - - τ τ)) \end{matrix}] \end{matrix} - - - - - - ((22)) . .$

优选地，在步骤2中，所述状态空间方程，如公式(3)所示：Preferably, in step 2, the state space equation, as shown in formula (3):

$\begin{matrix} [\begin{matrix} {\overset{\cdot \cdot}{z z}}_{11} ((t t)) \\ {\overset{\cdot \cdot}{z z}}_{22} ((t t)) \\ {\overset{\cdot\cdot \cdot\cdot}{z z}}_{11} ((t t)) \\ {\overset{\cdot\cdot \cdot\cdot}{z z}}_{22} ((t t)) \end{matrix}] = = [\begin{matrix} 00 & 00 & 11 & 00 \\ 00 & 00 & 00 & 11 \\ - - {ω ω}_{n no 11}^{22} - - {K K}_{f f} \frac{{ω ω}_{n no 11}^{22}}{{k k}_{11}} {a a}_{11} & 00 & - - 22 {ζ ζ}_{11} {ω ω}_{n no 11} & 00 \\ 00 & - - {ω ω}_{n no 22}^{22} - - {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} {a a}_{22} & 00 & - - 22 {ζ ζ}_{22} {ω ω}_{n no 22} \end{matrix}] [\begin{matrix} {z z}_{11} ((t t)) \\ {z z}_{22} ((t t)) \\ {\overset{\cdot \cdot}{z z}}_{11} ((t t)) \\ {\overset{\cdot &Center Dot;}{z z}}_{22} ((t t)) \end{matrix}] + + \\ [\begin{matrix} 00 & 00 & 00 & 00 \\ 00 & 00 & 00 & 00 \\ 00 & {K K}_{f f} \frac{{ω ω}_{n no 11}^{22}}{{k k}_{11}} {a a}_{11} & 00 & 00 \\ {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} {a a}_{11} & 00 & 00 & 00 \end{matrix}] [\begin{matrix} {z z}_{11} ((t t - - \frac{τ τ}{22})) \\ {z z}_{22} ((t t - - \frac{τ τ}{22})) \\ {\overset{\cdot &Center Dot;}{z z}}_{11} ((t t - - \frac{τ τ}{22})) \\ {\overset{\cdot \cdot}{z z}}_{22} ((t t - - \frac{τ τ}{22})) \end{matrix}] + + [\begin{matrix} 00 & 00 & 00 & 00 \\ 00 & 00 & 00 & 00 \\ 00 & 00 & 00 & 00 \\ 00 & {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} (({a a}_{22} - - {a a}_{11})) & 00 & 00 \end{matrix}] [\begin{matrix} {z z}_{11} ((t t - - τ τ)) \\ {z z}_{22} ((t t - - τ τ)) \\ {\overset{\cdot &Center Dot;}{z z}}_{11} ((t t - - τ τ)) \\ {\overset{\cdot &Center Dot;}{z z}}_{22} ((t t - - τ τ)) \end{matrix}] \end{matrix} - - - - - - ((33))$

令则公式(3)变为：make Then formula (3) becomes:

$\overset{\cdot &Center Dot;}{Z Z} ((t t)) = = A A Z Z ((t t)) + + {B B}_{11} Z Z ((t t - - \frac{τ τ}{22})) + + {B B}_{22} Z Z ((t t - - τ τ)) - - - - - - ((44))$

其中，表示状态速度，A表示状态位移项系数矩阵，B₁表示时滞状态位移项系数矩阵，B₂表示时滞状态位移项系数矩阵。in, Represents the state velocity, A represents the coefficient matrix of the state displacement term, B ₁ represents the coefficient matrix of the time-delay state displacement term, and B ₂ represents the coefficient matrix of the time-delay state displacement term.

优选地，所述步骤3，具体为：Preferably, the step 3 is specifically:

令t＝ξ·τ，则 Let t=ξ·τ, then

归一化后的状态空间方程变为标准形式，如公式(5)所示：The normalized state space equation becomes standard form, as shown in formula (5):

$\overset{\cdot &Center Dot;}{Z Z} ((ξ ξ)) = = τ τ A A Z Z ((ξ ξ)) + + {τB τB}_{11} Z Z ((ξ ξ - - 0.5 0.5)) + + {τB τB}_{22} Z Z ((ξ ξ - - 11)) - - - - - - ((55)) . .$

优选地，所述步骤4，具体为：Preferably, the step 4 is specifically:

在区间内取n+1个第二类切比雪夫离散点ξ_i，如公式(6)所示：Take n+1 second-type Chebyshev discrete points ξ _i in the interval, as shown in formula (6):

${ξ ξ}_{i i} = = \frac{11}{22} [[11 - - c c o o s the s ((\frac{i i π π}{n no}))]],, i i = = 00,, ... ...,, n no - - - - - - ((66))$

则在每一个离散点上都需满足方程(5)，即：Then equation (5) needs to be satisfied at each discrete point, namely:

$\overset{\cdot \cdot}{Z Z} (({ξ ξ}_{i i})) = = τ τ A A Z Z (({ξ ξ}_{i i})) + + {τB τB}_{11} Z Z (({ξ ξ}_{i i} - - 0.5 0.5)) + + {τB τB}_{22} Z Z (({ξ ξ}_{i i} - - 11)) - - - - - - ((77))$

则在[0,1]区间上，Then on the [0,1] interval,

$[\begin{matrix} \overset{\cdot &Center Dot;}{Z Z} (({ξ ξ}_{00})) \\ \overset{\cdot &Center Dot;}{Z Z} (({ξ ξ}_{11})) \\ . . \\ . . \\ . . \\ \overset{\cdot \cdot}{Z Z} (({ξ ξ}_{n no})) \end{matrix}] = = τ τ \cdot \cdot {I I}_{n no + + 11} &CircleTimes; &CircleTimes; A A [\begin{matrix} Z Z (({ξ ξ}_{00})) \\ Z Z (({ξ ξ}_{11})) \\ . . \\ . . \\ . . \\ Z Z (({ξ ξ}_{n no})) \end{matrix}] + + τ τ \cdot &Center Dot; {I I}_{n no + + 11} &CircleTimes; &CircleTimes; {B B}_{11} [\begin{matrix} Z Z (({ξ ξ}_{00} - - 0.5 0.5)) \\ Z Z (({ξ ξ}_{11} - - 0.5 0.5)) \\ . . \\ . . \\ . . \\ Z Z (({ξ ξ}_{n no} - - 0.5 0.5)) \end{matrix}] + + τ τ \cdot &Center Dot; {I I}_{n no + + 11} &CircleTimes; &CircleTimes; {B B}_{22} [\begin{matrix} Z Z (({ξ ξ}_{00} - - 11)) \\ Z Z (({ξ ξ}_{11} - - 11)) \\ . . \\ . . \\ . . \\ Z Z (({ξ ξ}_{n no} - - 11)) \end{matrix}] - - - - - - ((88))$

其中，表示Kronecker乘积，即I_n+1表示(n+1)维方阵。in, Denotes the Kronecker product, ie I _n+1 means (n+1) dimensional square matrix.

优选地，所述步骤5，具体为：Preferably, the step 5 is specifically:

利用微分求积法，用位移项表示公式(8)左端的速度项；为了表述方便，记为f(ξ_i)表示状态位移项标量表述符号；Using the differential quadrature method, the velocity term on the left side of formula (8) is represented by the displacement term; for the convenience of expression, record for f(ξ _i ) represents the scalar expression symbol of the state displacement item;

首先用(n+1)个点(ξ₀,f(ξ₀)),(ξ₁,f(ξ₁)),…,(ξ_n,f(ξ_n))进行拉格朗日插值，结果如公式(9)所示：First use (n+1) points (ξ ₀ ,f(ξ ₀ )),(ξ ₁ ,f(ξ ₁ )),…,(ξ _n ,f(ξ _n )) for Lagrangian interpolation, The result is shown in formula (9):

$f f ((ξ ξ)) = = {Σ Σ}_{j j = = 00}^{n no} {l l}_{j j} ((ξ ξ)) f f (({ξ ξ}_{j j})) - - - - - - ((99))$

其中，插值基函数为：Among them, the interpolation basis function is:

${l l}_{i i} ((ξ ξ)) = = \underset{k k &NotEqual; &NotEqual; i i}{{Π Π}_{k k = = 00}^{n no}} \frac{((ξ ξ - - {ξ ξ}_{k k}))}{(({ξ ξ}_{i i} - - {ξ ξ}_{k k}))} - - - - - - ((1010))$

对拉格朗日插值函数求导，再将公式(8)左端速度项各时间点代入，结果如公式(11)所示：Deriving the Lagrangian interpolation function, and then substituting the time points of the velocity item at the left end of formula (8), the result is shown in formula (11):

其中，H的表达式如公式(12)所示：Among them, the expression of H is shown in formula (12):

${H h}_{i i j j} = = \frac{{dl dl}_{j j} ((ξ ξ))}{d d ξ ξ} {| |}_{ξ ξ = = {ξ ξ}_{i i}} = = \{\begin{matrix} \frac{\underset{k k &NotEqual; &NotEqual; i i,, j j}{{Π Π}_{k k = = 00}^{n no}} (({ξ ξ}_{i i} - - {ξ ξ}_{k k}))}{\underset{k k &NotEqual; &NotEqual; j j}{{Π Π}_{k k = = 00}^{n no}} (({ξ ξ}_{j j} - - {ξ ξ}_{k k}))},, & j j &NotEqual; &NotEqual; i i \\ \underset{k k &NotEqual; &NotEqual; i i}{{Σ Σ}_{k k = = 00}^{n no}} \frac{11}{{ξ ξ}_{j j} - - {ξ ξ}_{k k}},, & j j = = i i \end{matrix} - - - - - - ((1212)) . .$

优选地，所述步骤6，具体为：Preferably, the step 6 is specifically:

时滞位移项Z(ξ-0.5)以[-1,0]和[0,1]区间上的第二类切比雪夫点表示；找到i使得ξ_i-0.5＜0且ξ_i+1-0.5＞0(i＝0,1,…,n)；以[-1,0]和[0,1]区间上的第二类切比雪夫点为插值点得到拉格朗日插值函数，再将时滞位移项的时间点代入，得到：The delay displacement term Z(ξ-0.5) is represented by a Chebyshev point of the second kind on the interval [-1,0] and [0,1]; find i such that ξ _i -0.5<0 and ξ _i+1 - 0.5>0 (i=0,1,...,n); use the second kind of Chebyshev points on the interval [-1,0] and [0,1] as the interpolation point to obtain the Lagrangian interpolation function, and then Substituting the time point of the time-delay displacement term, we get:

其中，in,

${T T}_{i i j j}^{d d} = = \{\begin{matrix} \frac{\underset{k k &NotEqual; &NotEqual; j j}{{Π Π}_{k k = = 00}^{n no}} [[(({ξ ξ}_{i i} - - 0.5 0.5)) - - (({ξ ξ}_{k k} - - 11))]]}{\underset{k k &NotEqual; &NotEqual; j j}{{Π Π}_{k k = = 00}^{n no}} [[(({ξ ξ}_{j j} - - 11)) - - (({ξ ξ}_{k k} - - 11))]]},, & (({ξ ξ}_{i i} - - 0.5 0.5)) &NotElement; &NotElement; {{{ξ ξ}_{k k} - - 11 | | k k = = 00,, 11,, ... ...,, n no}} \\ 00,, & (({ξ ξ}_{i i} - - 0.5 0.5)) &Element; &Element; {{{ξ ξ}_{k k} - - 11 | | k k = = 00,, 11,, ... ...,, j j - - 11,, j j + + 11,, ... ...,, n no}} \\ 11,, & {ξ ξ}_{i i} - - 0.5 0.5 = = {ξ ξ}_{j j} - - 11 \end{matrix} - - - - - - ((1414))$

其中，in,

${T T}_{i i j j}^{p p} = = \{\begin{matrix} \frac{\underset{k k &NotEqual; &NotEqual; j j}{{Π Π}_{k k = = 00}^{n no}} [[(({ξ ξ}_{i i} - - 0.5 0.5)) - - {ξ ξ}_{k k}]]}{\underset{k k &NotEqual; &NotEqual; j j}{{Π Π}_{k k = = 00}^{n no}} (({ξ ξ}_{i i} - - {ξ ξ}_{k k}))},, & (({ξ ξ}_{i i} - - 0.5 0.5)) &NotElement; &NotElement; {{{ξ ξ}_{k k} - - 11 | | k k = = 00,, 11,, ... ...,, n no}} \\ 00,, & (({ξ ξ}_{i i} - - {0.5 0.5}_{11})) &Element; &Element; {{{ξ ξ}_{k k} | | k k = = 00,, 11,, ... ...,, j j - - 11,, j j + + 11,, ... ...,, n no}} \\ 11,, & {ξ ξ}_{i i} - - 0.5 0.5 = = {ξ ξ}_{j j} \end{matrix} - - - - - - ((1616))$

时滞项Z(ξ-1)在区间[-1,0]上的离散，如公式(17)所示：The discretization of the delay term Z(ξ-1) on the interval [-1,0] is shown in formula (17):

能够得到，T＝I_n+1，It can be obtained that T=I _n+1 ,

其中，T表示 Among them, T means

优选地，所述步骤7，具体为：Preferably, the step 7 is specifically:

构造Floquet传递矩阵：Construct the Floquet transfer matrix:

矩阵H,T^d,T^p,T消掉第一行后分别记为 The matrices H, T ^d , T ^p , T are recorded as

其中，T^d表示T^p表示 Among them, T ^d represents T ^p said

综合以上各式，可以得出相邻两个区间[-1,0]和[0,1]之间的状态转移关系为：Combining the above formulas, it can be concluded that the state transition relationship between two adjacent intervals [-1,0] and [0,1] is:

$P P [\begin{matrix} Z Z (({ξ ξ}_{00})) \\ Z Z (({ξ ξ}_{11})) \\ . . \\ . . \\ . . \\ Z Z (({ξ ξ}_{n no})) \end{matrix}] = = Q Q [\begin{matrix} Z Z (({ξ ξ}_{00} - - 11)) \\ Z Z (({ξ ξ}_{11} - - 11)) \\ . . \\ . . \\ . . \\ Z Z (({ξ ξ}_{n no} - - 11)) \end{matrix}] - - - - - - ((1818))$

其中，in,

Floquet状态转移矩阵Φ如公式所示：The Floquet state transition matrix Φ is shown in the formula:

其中，表示Penrose-Moor广义逆；in, Represents the Penrose-Moor generalized inverse;

根据Floquet理论，若Φ的所有特征值的模均小于1，则系统是稳定的；若Φ中任一特征值的模大于1，则系统是不稳定的。According to Floquet theory, if the modulus of all eigenvalues of Φ are less than 1, the system is stable; if the modulus of any eigenvalue of Φ is greater than 1, the system is unstable.

优选地，还包括如下步骤：Preferably, the following steps are also included:

步骤8：并绘制出系统在时滞参数空间的稳定性图谱。Step 8: And draw the stability map of the system in the delay parameter space.

与现有技术相比，本发明具有如下的有益效果：Compared with the prior art, the present invention has the following beneficial effects:

本发明与传统单车刀车削加工相比，采用微分求积法分析多车刀并行车削系统动力学特性，获得优化后的切削参数，极大地提高了加工效率。Compared with the traditional single-tool turning process, the present invention adopts the differential quadrature method to analyze the dynamic characteristics of the multi-tool parallel turning system, obtains optimized cutting parameters, and greatly improves the machining efficiency.

附图说明Description of drawings

通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：Other characteristics, objects and advantages of the present invention will become more apparent by reading the detailed description of non-limiting embodiments made with reference to the following drawings:

图1为多车刀并行车削系统示意图，两把车刀安装在不同的刀架上。Figure 1 is a schematic diagram of a parallel turning system with multiple turning tools. Two turning tools are installed on different tool holders.

图2为多车刀并行车削在车刀切削深度参数空间(a₁,a₂)的稳定性图谱，其中黑色部分表示不稳定区域。Fig. 2 is the stability map of parallel turning with multiple turning tools in the cutting depth parameter space (a ₁ , a ₂ ), where the black part represents the unstable region.

图1中，各代式所表示的含义如下：In Figure 1, the meanings represented by each generation formula are as follows:

A表示第一车刀，B表示第二车刀，C_w表示工件的阻尼系数，K_w表示工件的刚度系数，C₁表示第一车刀的阻尼系数，K₁表示第一车刀的刚度系数，C₂表示第二车刀的阻尼系数，K₂表示第二车刀的刚度系数，a₁表示第一车刀的切削深度，a₂表示第二车刀的切削深度，z₁表示第一车刀的动力学响应位移，z₂表示第二车刀的动力学响应位移，Z表示车刀柔度方向。A represents the first turning tool, B represents the second turning tool, C _w represents the damping coefficient of the workpiece, K _w represents the stiffness coefficient of the workpiece, C ₁ represents the damping coefficient of the first turning tool, K ₁ represents the stiffness of the first turning tool coefficient, C ₂ represents the damping coefficient of the second turning tool, K ₂ represents the stiffness coefficient of the second turning tool, a ₁ represents the cutting depth of the first turning tool, a ₂ represents the cutting depth of the second turning tool, z ₁ represents the cutting depth of the second turning tool The dynamic response displacement of the first turning tool, z ₂ represents the dynamic response displacement of the second turning tool, and Z represents the flexibility direction of the turning tool.

图2中，a₁表示第一车刀的切削深度，a₂表示第二车刀的切削深度。In Fig. 2, a ₁ represents the cutting depth of the first turning tool, and a ₂ represents the cutting depth of the second turning tool.

具体实施方式detailed description

下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进。这些都属于本发明的保护范围。The present invention will be described in detail below in conjunction with specific embodiments. The following examples will help those skilled in the art to further understand the present invention, but do not limit the present invention in any form. It should be noted that those skilled in the art can make several modifications and improvements without departing from the concept of the present invention. These all belong to the protection scope of the present invention.

本发明提供的基于微分求积法的多车刀并行车削稳定性判定方法，包括以下步骤：The multi-turn tool parallel turning stability determination method based on the differential quadrature method provided by the present invention comprises the following steps:

(1)建立多车刀并行车削加工的动力学方程，即多时滞二阶微分方程。根据振动试验/模态试验结果确定方程中的相关参数。(1) Establish the dynamic equation of multi-tool parallel turning, that is, the multi-time-delay second-order differential equation. Determine the relevant parameters in the equation based on the vibration test/modal test results.

(2)对微分方程进行状态空间变换，得到状态空间方程。(2) Perform state space transformation on the differential equation to obtain the state space equation.

(3)对状态空间方程进行归一化处理，将其化为便于微分求积法求解的标准形式。(3) Normalize the state-space equations and turn them into a standard form that is easy to solve by the differential quadrature method.

(4)基于微分求积法，用相邻单位区间[0,1]或者[-1,0]上的离散点表示状态空间方程中的导数项及时滞项，将多时滞方程等价转化为一组只含切削参数的代数方程组。(4) Based on the differential quadrature method, the discrete points on the adjacent unit interval [0,1] or [-1,0] are used to represent the derivative term and the delay term in the state space equation, and the multi-delay equation is equivalently transformed into A set of algebraic equations containing only cutting parameters.

(5)基于Floquet定理判定多车刀并行车削加工过程稳定性，并在切削参数空间画出稳定性图谱。(5) Based on the Floquet theorem, the stability of the multi-tool parallel turning process is determined, and the stability map is drawn in the cutting parameter space.

更为具体地，下面结合具体加工实例说明本发明的具体实施方案，实例参数引自文献2中的实验2.采用两把车刀并行车削，且车刀安装于不同的刀架上，相互之间没有振动耦合，其示意图如图1所示。记切深较小的车刀为第一车刀，切深较大的车刀为第二车刀，加工工件的直径d为d＝35mm，主轴转速Ω为Ω＝1800rpm，进给方向切削力系数K_f为K_f＝872MPa，模态试验参数如表1所示：More specifically, the specific implementation of the present invention will be described below in conjunction with specific processing examples. The example parameters are quoted from Experiment 2 in Document 2. Two turning tools are used for parallel turning, and the turning tools are installed on different tool holders. There is no vibration coupling between them, as shown in Figure 1. Note that the turning tool with smaller cutting depth is the first turning tool, and the turning tool with larger cutting depth is the second turning tool. The diameter d of the workpiece to be processed is d=35mm, the spindle speed Ω is Ω=1800rpm, and the cutting force in the feed direction The coefficient K _f is K _f =872MPa, and the modal test parameters are shown in Table 1:

表1多车刀并行切削实例模态试验参数Table 1 Modal test parameters of multi-tool parallel cutting example

模态振型Mode shape f_n(Hz)f _n (Hz) ζ(％)ζ(%) k(N/m)k(N/m) 11 2238.92238.9 3.233.23 4.769*10⁷ 4.769*10 ⁷ 22 2372.32372.3 4.514.51 1.166*10⁸ 1.166*10 ⁸

其中，f_n表示固有频率，ζ表示相对阻尼比，k表示刚度系数。Among them, f _n represents the natural frequency, ζ represents the relative damping ratio, and k represents the stiffness coefficient.

步骤(1)，建立多车刀并行车削动力学方程，如公式(1)所示：In step (1), establish the multi-tool parallel turning dynamic equation, as shown in formula (1):

$\begin{matrix} [\begin{matrix} {\overset{\cdot\cdot \cdot\cdot}{z z}}_{11} ((t t)) \\ {\overset{\cdot\cdot \cdot\cdot}{z z}}_{22} ((t t)) \end{matrix}] + + [\begin{matrix} 22 {ζ ζ}_{11} {ω ω}_{n no 11} & 00 \\ 00 & 22 {ζ ζ}_{22} {ω ω}_{n no 22} \end{matrix}] [\begin{matrix} {\overset{\cdot &Center Dot;}{z z}}_{11} ((t t)) \\ {\overset{\cdot &Center Dot;}{z z}}_{22} ((t t)) \end{matrix}] + + [\begin{matrix} {ω ω}_{n no 11}^{22} & 00 \\ 00 & {ω ω}_{n no 22}^{22} \end{matrix}] [\begin{matrix} {z z}_{11} ((t t)) \\ {z z}_{22} ((t t)) \end{matrix}] = = \\ {K K}_{f f} [\begin{matrix} \frac{{ω ω}_{n no 11}^{22}}{{k k}_{11}} {a a}_{11} ((- - {z z}_{11} ((t t)) + + {z z}_{22} ((t t - - \frac{τ τ}{22})))) \\ \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} {a a}_{11} ((- - {z z}_{22} ((t t)) + + {z z}_{11} ((t t - - \frac{τ τ}{22})))) + + \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} (({a a}_{22} - - {a a}_{11})) ((- - {z z}_{22} ((t t)) + + {z z}_{22} ((t t - - τ τ)))) \end{matrix}] \end{matrix} - - - - - - ((11))$

其中，下标1表示第一车刀，下标2表示第二车刀，z(t)表示动力学响应位移，ω_n表示固有圆频率，K_f表示进给方向的切削力系数，τ表示工件旋转周期，ζ表示相对阻尼比，k表示刚度系数，a₁表示第一车刀的切削深度，a₂表示第二车刀的切削深度，表示加速度，表示速度,t表示时间。Among them, subscript 1 indicates the first turning tool, subscript 2 indicates the second turning tool, z(t) indicates the dynamic response displacement, ω _n indicates the natural circular frequency, K _f indicates the cutting force coefficient in the feed direction, and τ indicates The workpiece rotation period, ζ represents the relative damping ratio, k represents the stiffness coefficient, a ₁ represents the cutting depth of the first turning tool, a ₂ represents the cutting depth of the second turning tool, represents the acceleration, means speed, and t means time.

对公式(1)进行整理，得多时滞二阶微分方程如公式(2)所示：Arranging formula (1), the second-order differential equation with many delays is shown in formula (2):

$\begin{matrix} [\begin{matrix} {\overset{\cdot\cdot \cdot\cdot}{z z}}_{11} ((t t)) \\ {\overset{\cdot\cdot \cdot\cdot}{z z}}_{22} ((t t)) \end{matrix}] + + [\begin{matrix} 22 {ζ ζ}_{11} {ω ω}_{11 n no} & 00 \\ 00 & 22 {ζ ζ}_{22} {ω ω}_{n no 22} \end{matrix}] [\begin{matrix} {\overset{\cdot &Center Dot;}{z z}}_{11} ((t t)) \\ {\overset{\cdot &Center Dot;}{z z}}_{22} ((t t)) \end{matrix}] + + [\begin{matrix} {ω ω}_{n no 11}^{22} + + {K K}_{f f} \frac{{ω ω}_{n no 11}^{22}}{{k k}_{11}} {a a}_{11} & 00 \\ 00 & {ω ω}_{n no 22}^{22} + + {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} {a a}_{22} \end{matrix}] [\begin{matrix} {z z}_{11} ((t t)) \\ {z z}_{22} ((t t)) \end{matrix}] = = \\ [\begin{matrix} 00 & {K K}_{f f} \frac{{ω ω}_{n no 11}^{22}}{{k k}_{11}} {a a}_{11} \\ {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} {a a}_{11} & 00 \end{matrix}] [\begin{matrix} {z z}_{11} ((t t - - \frac{τ τ}{22})) \\ {z z}_{22} ((t t - - \frac{τ τ}{22})) \end{matrix}] + + [\begin{matrix} 00 & 00 \\ 00 & {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} (({a a}_{22} - - {a a}_{11})) \end{matrix}] [\begin{matrix} {z z}_{11} ((t t - - τ τ)) \\ {z z}_{22} ((t t - - τ τ)) \end{matrix}] \end{matrix} - - - - - - ((22))$

步骤(2)，对公式(2)进行状态空间变换，得状态空间方程如公式(3)所示：In step (2), the state space transformation is performed on formula (2), and the state space equation is obtained as shown in formula (3):

$\begin{matrix} [\begin{matrix} {\overset{\cdot \cdot}{z z}}_{11} ((t t)) \\ {\overset{\cdot \cdot}{z z}}_{22} ((t t)) \\ {\overset{\cdot\cdot \cdot\cdot}{z z}}_{11} ((t t)) \\ {\overset{\cdot\cdot \cdot\cdot}{z z}}_{22} ((t t)) \end{matrix}] = = [\begin{matrix} 00 & 00 & 11 & 00 \\ 00 & 00 & 00 & 11 \\ - - {ω ω}_{n no 11}^{22} - - {K K}_{f f} \frac{{ω ω}_{n no 11}^{22}}{{k k}_{11}} {a a}_{11} & 00 & - - 22 {ζ ζ}_{11} {ω ω}_{n no 11} & 00 \\ 00 & - - {ω ω}_{n no 22}^{22} - - {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} {a a}_{22} & 00 & - - 22 {ζ ζ}_{22} {ω ω}_{n no 22} \end{matrix}] [\begin{matrix} {z z}_{11} ((t t)) \\ {z z}_{22} ((t t)) \\ {\overset{\cdot \cdot}{z z}}_{11} ((t t)) \\ {\overset{\cdot \cdot}{z z}}_{22} ((t t)) \end{matrix}] + + \\ [\begin{matrix} 00 & 00 & 00 & 00 \\ 00 & 00 & 00 & 00 \\ 00 & {K K}_{f f} \frac{{ω ω}_{n no 11}^{22}}{{k k}_{11}} {a a}_{11} & 00 & 00 \\ {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} {a a}_{11} & 00 & 00 & 00 \end{matrix}] [\begin{matrix} {z z}_{11} ((t t - - \frac{τ τ}{22})) \\ {z z}_{22} ((t t - - \frac{τ τ}{22})) \\ {\overset{\cdot &Center Dot;}{z z}}_{11} ((t t - - \frac{τ τ}{22})) \\ {\overset{\cdot &Center Dot;}{z z}}_{22} ((t t - - \frac{τ τ}{22})) \end{matrix}] + + [\begin{matrix} 00 & 00 & 00 & 00 \\ 00 & 00 & 00 & 00 \\ 00 & 00 & 00 & 00 \\ 00 & {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} (({a a}_{22} - - {a a}_{11})) & 00 & 00 \end{matrix}] [\begin{matrix} {z z}_{11} ((t t - - τ τ)) \\ {z z}_{22} ((t t - - τ τ)) \\ {\overset{\cdot &Center Dot;}{z z}}_{11} ((t t - - τ τ)) \\ {\overset{\cdot &Center Dot;}{z z}}_{22} ((t t - - τ τ)) \end{matrix}] \end{matrix} - - - - - - ((33))$

令则公式(3)变为：make Then formula (3) becomes:

其中，表示状态速度项，A表示状态位移项系数矩阵，B₁表示时滞状态位移项系数矩阵，B₂表示时滞状态位移项系数矩阵。in, Represents the state velocity item, A represents the coefficient matrix of the state displacement term, B ₁ represents the coefficient matrix of the time-delay state displacement term, and B ₂ represents the coefficient matrix of the time-delay state displacement term.

步骤(3)，对状态空间方程进行归一化处理。不妨令t＝ξ·τ，则归一化后的状态空间方程变为标准形式，如公式(5)所示：Step (3), normalize the state space equation. Let t=ξ·τ, then The normalized state space equation becomes standard form, as shown in formula (5):

$\overset{\cdot \cdot}{Z Z} ((ξ ξ)) = = τ τ A A Z Z ((ξ ξ)) + + {τB τB}_{11} Z Z ((ξ ξ - - 0.5 0.5)) + + {τB τB}_{22} Z Z ((ξ ξ - - 11)) - - - - - - ((55))$

步骤(4)，对公式(5)进行周期离散，将其等价转化为一组代数方程。在区间内取n+1个第二类切比雪夫离散点ξ_i，如公式(6)所示：In step (4), perform periodic discretization on the formula (5), and transform it into a set of algebraic equations equivalently. Take n+1 second-type Chebyshev discrete points ξ _i in the interval, as shown in formula (6):

则在[0,1]区间上，Then on the [0,1] interval,

$[\begin{matrix} \overset{\cdot \cdot}{Z Z} (({ξ ξ}_{00})) \\ \overset{\cdot \cdot}{Z Z} (({ξ ξ}_{11})) \\ . . \\ . . \\ . . \\ \overset{\cdot \cdot}{Z Z} (({ξ ξ}_{n no})) \end{matrix}] = = τ τ \cdot \cdot {I I}_{n no + + 11} &CircleTimes; &CircleTimes; A A [\begin{matrix} Z Z (({ξ ξ}_{00})) \\ Z Z (({ξ ξ}_{11})) \\ . . \\ . . \\ . . \\ Z Z (({ξ ξ}_{n no})) \end{matrix}] + + τ τ \cdot \cdot {I I}_{n no + + 11} &CircleTimes; &CircleTimes; {B B}_{11} [\begin{matrix} Z Z (({ξ ξ}_{00} - - 0.5 0.5)) \\ Z Z (({ξ ξ}_{11} - - 0.5 0.5)) \\ . . \\ . . \\ . . \\ Z Z (({ξ ξ}_{n no} - - 0.5 0.5)) \end{matrix}] + + τ τ \cdot &Center Dot; {I I}_{n no + + 11} &CircleTimes; &CircleTimes; {B B}_{22} [\begin{matrix} Z Z (({ξ ξ}_{00} - - 11)) \\ Z Z (({ξ ξ}_{11} - - 11)) \\ . . \\ . . \\ . . \\ Z Z (({ξ ξ}_{n no} - - 11)) \end{matrix}] - - - - - - ((88))$

步骤(5)，微分求积法Step (5), differential quadrature method

利用微分求积法，用位移项表示公式(8)左端的速度项。为了表述方便，记为f(ξ_i)表示状态位移项标量表述符号。首先用(n+1)个点(ξ₀,f(ξ₀)),(ξ₁,f(ξ₁)),…,(ξ_n,f(ξ_n))进行拉格朗日(Lagrange)插值，结果如公式(9)所示：Using the differential quadrature method, the velocity term at the left end of formula (8) is represented by the displacement term. For the convenience of expression, remember for f(ξ _i ) represents the scalar expression symbol of the state displacement item. First use (n+1) points (ξ ₀ , f(ξ ₀ )), (ξ ₁ , f(ξ ₁ )),…, (ξ _n , f(ξ _n )) to perform Lagrange ) interpolation, the result is shown in formula (9):

其中，插值基函数为：Among them, the interpolation basis function is:

${H h}_{i i j j} = = \frac{{dl dl}_{j j} ((ξ ξ))}{d d ξ ξ} {| |}_{ξ ξ = = {ξ ξ}_{i i}} = = \{\begin{matrix} \frac{\underset{k k &NotEqual; &NotEqual; i i,, j j}{{Π Π}_{k k = = 00}^{n no}} (({ξ ξ}_{i i} - - {ξ ξ}_{k k}))}{\underset{k k &NotEqual; &NotEqual; j j}{{Π Π}_{k k = = 00}^{n no}} (({ξ ξ}_{j j} - - {ξ ξ}_{k k}))},, & j j &NotEqual; &NotEqual; i i \\ \underset{k k &NotEqual; &NotEqual; i i}{{Σ Σ}_{k k = = 00}^{n no}} \frac{11}{{ξ ξ}_{j j} - - {ξ ξ}_{k k}},, & j j = = i i \end{matrix} - - - - - - ((1212))$

时滞位移项Z(ξ-0.5)要以[-1,0]和[0,1]区间上的第二类切比雪夫点表示。可以找到i使得ξ_i-0.5＜0且ξ_i+1-0.5＞0(i＝0,1,…,n)；以[-1,0]和[0,1]区间上的第二类切比雪夫点为插值点得到拉格朗日插值函数，再将时滞位移项的时间点代入，可得：The time-delay displacement term Z(ξ-0.5) should be represented by Chebyshev points of the second kind on the interval [-1,0] and [0,1]. i can be found such that ξ _i -0.5<0 and ξ _i+1 -0.5>0 (i=0,1,...,n); the second category on the interval [-1,0] and [0,1] The Chebyshev point is used as the interpolation point to obtain the Lagrangian interpolation function, and then the time point of the time-delay displacement item is substituted to obtain:

其中，in,

时滞项Z(ξ-1)在区间[-1,0]上的离散十分简便，如公式(17)所示：The discretization of the delay term Z(ξ-1) on the interval [-1,0] is very simple, as shown in formula (17):

易知，T＝I_n+1，It is easy to know that T=I _n+1 ,

其中，T表示 Among them, T means

步骤(6)，构造Floquet传递矩阵Step (6), constructing the Floquet transfer matrix

其中，T^d表示T^p表示 Among them, T ^d represents T ^p said

其中，in,

在第一车刀和第二车刀的切削深度参数空间(a₁,a₂)，多车刀并行车削的稳定性图谱如图2所示。In the cutting depth parameter space (a ₁ , a ₂ ) of the first turning tool and the second turning tool, the stability map of parallel turning with multiple turning tools is shown in Fig. 2 .

根据文献1，若保持其他实验条件不变，只采用第一车刀进行传统车削，稳定性加工的极限切深为3.4mm，只采用第二车刀进行传统车削，稳定性加工的极限切深为12.6mm。根据附图2，采用多车刀并行车削之后，稳定性加工的极限切深获得极大提高，若采用(a₁,a₂)＝(8,26)切深参数组合进行无再生颤振稳定性加工，第一车刀的切深由3.4mm提高到8mm，提高了135.3％，第二车刀的切深由12.6mm提高到26mm，提高了106.3％。由此可见，与传统单车刀切削相比，多车刀并行车削极大地提高了加工效率，加工效率的提高不仅体现在多把车刀同时车削上，还体现在多车刀并行车削加工极大地提高了传统单车刀车削加工的稳定区域。According to Document 1, if other experimental conditions remain unchanged, only the first turning tool is used for traditional turning, the limit depth of cut for stable machining is 3.4mm, and only the second turning tool is used for traditional turning, the limit depth of cut for stable machining It is 12.6mm. According to Figure 2, after parallel turning with multiple turning tools, the limit depth of cut for stable machining is greatly improved. If (a ₁ ,a ₂ )=(8,26) depth of cut For permanent processing, the cutting depth of the first turning tool is increased from 3.4mm to 8mm, which is an increase of 135.3%, and the cutting depth of the second turning tool is increased from 12.6mm to 26mm, which is an increase of 106.3%. It can be seen that compared with the traditional single turning tool cutting, the parallel turning of multiple turning tools greatly improves the processing efficiency. Increased stable area for conventional single-tool turning operations.

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变形或修改，这并不影响本发明的实质内容。Specific embodiments of the present invention have been described above. It should be understood that the present invention is not limited to the specific embodiments described above, and those skilled in the art may make various changes or modifications within the scope of the claims, which do not affect the essence of the present invention.

Claims

1. a method for judging the stability of parallel turning of many turning tools based on differential quadrature method, is characterized in that, comprises the steps:

Step 1: Establish the multi-tool parallel turning dynamic equation, sort out the multi-turn tool parallel turning dynamic equation, and multi-time-delay second-order differential equation;

Step 2: performing state-space transformation on the multi-delay second-order differential equation to obtain a state-space equation;

Step 3: performing normalization processing on the state-space equation to obtain the state-space equation in standard form;

Step 4: Periodically discretize the state-space equation in standard form, and convert it into a set of algebraic equations as the expression of the state-space equation;

Step 5: Based on the Lagrangian interpolation function, take the second kind of Chebyshev point on the unit interval [0,1] as a discrete point, and use the differential quadrature method to represent the derivative in the expression of the state space equation with the displacement term item;

Step 6: Based on the Lagrangian interpolation function, judge the interval of the delay item in the expression of the state space equation. If the discrete point of the delay item belongs to the interval [0,1], then use [0,1] The second kind of Chebyshev point is a discrete point to represent the delay term; if the discrete point of the delay term belongs to the interval [-1,0], then the second kind of Chebyshev point on [-1,0] is Discrete points represent delay items;

Step 7: Construct the state transition matrix between adjacent unit intervals [0,1] and [-1,0], and judge the stability of the original system according to the Floquet theory.

2. the method for judging the stability of multi-turning tool parallel turning based on differential quadrature method according to claim 1, is characterized in that, in described step 1, described multi-turning tool parallel turning kinetic equation, as formula (1 ) as shown:

\begin{matrix} [\begin{matrix} {\overset{\cdot\cdot \cdot\cdot}{z z}}_{11} ((t t)) \\ {\overset{\cdot\cdot \cdot\cdot}{z z}}_{22} ((t t)) \end{matrix}] + + [\begin{matrix} 22 {ζ ζ}_{11} {ω ω}_{n no 11} & 00 \\ 00 & 22 {ζ ζ}_{22} {ω ω}_{n no 22} \end{matrix}] [\begin{matrix} {\overset{\cdot &Center Dot;}{z z}}_{11} ((t t)) \\ {\overset{\cdot &Center Dot;}{z z}}_{22} ((t t)) \end{matrix}] + + [\begin{matrix} {ω ω}_{n no 11}^{22} & 00 \\ 00 & {ω ω}_{n no 22}^{22} \end{matrix}] [\begin{matrix} {z z}_{11} ((t t)) \\ {z z}_{22} ((t t)) \end{matrix}] = = \\ {K K}_{f f} [\begin{matrix} \frac{{ω ω}_{n no 11}^{22}}{{k k}_{11}} {a a}_{11} ((- - {z z}_{11} ((t t)) + + {z z}_{22} ((t t - - \frac{τ τ}{22})))) \\ \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} {a a}_{11} ((- - {z z}_{22} ((t t)) + + {z z}_{11} ((t t - - \frac{τ τ}{22})))) + + \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} (({a a}_{22} - - {a a}_{11})) ((- - {z z}_{22} ((t t)) + + {z z}_{22} ((t t - - τ τ)))) \end{matrix}] \end{matrix} - - - - - - ((11))

Among them, subscript 1 indicates the first turning tool, subscript 2 indicates the second turning tool, z(t) indicates the dynamic response displacement, ω _n indicates the natural circular frequency, K _f indicates the cutting force coefficient in the feed direction, and τ indicates The workpiece rotation period, ζ represents the relative damping ratio, k represents the stiffness coefficient, a ₁ represents the cutting depth of the first turning tool, a ₂ represents the cutting depth of the second turning tool, represents the dynamic acceleration response, Indicates dynamic velocity response, t indicates time;

The multi-time-delay second-order differential equation, as shown in formula (2):

\begin{matrix} [\begin{matrix} {\overset{\cdot\cdot \cdot\cdot}{z z}}_{11} ((t t)) \\ {\overset{\cdot\cdot \cdot\cdot}{z z}}_{22} ((t t)) \end{matrix}] + + [\begin{matrix} 22 {ζ ζ}_{11} {ω ω}_{11 n no} & 00 \\ 00 & 22 {ζ ζ}_{22} {ω ω}_{n no 22} \end{matrix}] [\begin{matrix} {\overset{\cdot \cdot}{z z}}_{11} ((t t)) \\ {\overset{\cdot \cdot}{z z}}_{22} ((t t)) \end{matrix}] + + [\begin{matrix} {ω ω}_{n no 11}^{22} + + {K K}_{f f} \frac{{ω ω}_{n no 11}^{22}}{{k k}_{11}} {a a}_{11} & 00 \\ 00 & {ω ω}_{n no 22}^{22} + + {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} {a a}_{22} \end{matrix}] [\begin{matrix} {z z}_{11} ((t t)) \\ {z z}_{22} ((t t)) \end{matrix}] = = \\ [\begin{matrix} 00 & {K K}_{f f} \frac{{ω ω}_{n no 11}^{22}}{{k k}_{11}} {a a}_{11} \\ {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} {a a}_{11} & 00 \end{matrix}] [\begin{matrix} {z z}_{11} ((t t - - \frac{τ τ}{22})) \\ {z z}_{22} ((t t - - \frac{τ τ}{22})) \end{matrix}] + + [\begin{matrix} 00 & 00 \\ 00 & {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} (({a a}_{22} - - {a a}_{11})) \end{matrix}] [\begin{matrix} {z z}_{11} ((t t - - τ τ)) \\ {z z}_{22} ((t t - - τ τ)) \end{matrix}] \end{matrix} - - - - - - ((22)) . .

3. the method for judging stability of many turning tools parallel turning based on differential quadrature method according to claim 2, is characterized in that, in step 2, described state space equation, as shown in formula (3):

\begin{matrix} [\begin{matrix} {\overset{\cdot &Center Dot;}{z z}}_{11} ((t t)) \\ {\overset{\cdot &Center Dot;}{z z}}_{22} ((t t)) \\ {\overset{\cdot\cdot \cdot\cdot}{z z}}_{11} ((t t)) \\ {\overset{\cdot\cdot \cdot\cdot}{z z}}_{22} ((t t)) \end{matrix}] = = [\begin{matrix} 00 & 00 & 11 & 00 \\ 00 & 00 & 00 & 11 \\ - - {ω ω}_{n no 11}^{22} - - {K K}_{f f} \frac{{ω ω}_{n no 11}^{22}}{{k k}_{11}} {a a}_{11} & 00 & - - 22 {ζ ζ}_{11} {ω ω}_{n no 11} & 00 \\ 00 & - - {ω ω}_{n no 22}^{22} - - {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} {a a}_{22} & 00 & - - 22 {ζ ζ}_{22} {ω ω}_{n no 22} \end{matrix}] [\begin{matrix} {z z}_{11} ((t t)) \\ {z z}_{22} ((t t)) \\ {\overset{\cdot &Center Dot;}{z z}}_{11} ((t t)) \\ {\overset{\cdot &Center Dot;}{z z}}_{22} ((t t)) \end{matrix}] + + \\ [\begin{matrix} 00 & 00 & 00 & 00 \\ 00 & 00 & 00 & 00 \\ 00 & {K K}_{f f} \frac{{ω ω}_{n no 11}^{22}}{{k k}_{11}} {a a}_{11} & 00 & 00 \\ {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} {a a}_{11} & 00 & 00 & 00 \end{matrix}] [\begin{matrix} {z z}_{11} ((t t - - \frac{τ τ}{22})) \\ {z z}_{22} ((t t - - \frac{τ τ}{22})) \\ {\overset{\cdot &Center Dot;}{z z}}_{11} ((t t - - \frac{τ τ}{22})) \\ {\overset{\cdot &Center Dot;}{z z}}_{22} ((t t - - \frac{τ τ}{22})) \end{matrix}] + + [\begin{matrix} 00 & 00 & 00 & 00 \\ 00 & 00 & 00 & 00 \\ 00 & 00 & 00 & 00 \\ 00 & {K K}_{f f} \frac{{ω ω}_{n no 22}^{22}}{{k k}_{22}} (({a a}_{22} - - {a a}_{11})) & 00 & 00 \end{matrix}] [\begin{matrix} {z z}_{11} ((t t - - τ τ)) \\ {z z}_{22} ((t t - - τ τ)) \\ {\overset{\cdot \cdot}{z z}}_{11} ((t t - - τ τ)) \\ {\overset{\cdot \cdot}{z z}}_{22} ((t t - - τ τ)) \end{matrix}] \end{matrix} - - - - - - ((33))

make Then formula (3) becomes:

\overset{\cdot &Center Dot;}{Z Z} ((t t)) = = A A Z Z ((t t)) + + {B B}_{11} Z Z ((t t - - \frac{τ τ}{22})) + + {B B}_{22} Z Z ((t t - - τ τ)) - - - - - - ((44))

in, Represents the state velocity, A represents the coefficient matrix of the state displacement term, B ₁ represents the coefficient matrix of the time-delay state displacement term, and B ₂ represents the coefficient matrix of the time-delay state displacement term.

4. the multi-turn tool parallel turning stability determination method based on differential quadrature method according to claim 3, is characterized in that, described step 3, is specifically:

Let t=ξ·τ, then

The normalized state space equation becomes standard form, as shown in formula (5):

\overset{\cdot &Center Dot;}{Z Z} ((ξ ξ)) = = τ τ A A Z Z ((ξ ξ)) + + {τB τB}_{11} Z Z ((ξ ξ - - 0.5 0.5)) + + {τB τB}_{22} Z Z ((ξ ξ - - 11)) - - - - - - ((55)) . .

5. the multi-turn tool parallel turning stability determination method based on differential quadrature method according to claim 4, is characterized in that, described step 4, is specifically:

Take n+1 second-type Chebyshev discrete points ξ _i in the interval, as shown in formula (6):

{ξ ξ}_{i i} = = \frac{11}{22} [[11 - - c c o o s the s ((\frac{i i π π}{n no}))]],, i i = = 00,, 11,, ... ...,, n no - - - - - - ((66))

Then equation (5) needs to be satisfied at each discrete point, namely:

\overset{\cdot &Center Dot;}{Z Z} (({ξ ξ}_{i i})) = = τ τ A A Z Z (({ξ ξ}_{i i})) + + {τB τB}_{11} Z Z (({ξ ξ}_{i i} - - 0.5 0.5)) + + {τB τB}_{22} Z Z (({ξ ξ}_{i i} - - 11)) - - - - - - ((77))

Then on the [0,1] interval,

[\begin{matrix} \overset{\cdot \cdot}{Z Z} (({ξ ξ}_{00})) \\ \overset{\cdot \cdot}{Z Z} (({ξ ξ}_{11})) \\ . . \\ . . \\ . . \\ \overset{\cdot \cdot}{Z Z} (({ξ ξ}_{n no})) \end{matrix}] = = τ τ \cdot &Center Dot; {I I}_{n no + + 11} &CircleTimes; &CircleTimes; A A [\begin{matrix} Z Z (({ξ ξ}_{00})) \\ Z Z (({ξ ξ}_{11})) \\ . . \\ . . \\ . . \\ Z Z (({ξ ξ}_{n no})) \end{matrix}] + + τ τ \cdot \cdot {I I}_{n no + + 11} &CircleTimes; &CircleTimes; {B B}_{11} [\begin{matrix} Z Z (({ξ ξ}_{00} - - 0.5 0.5)) \\ Z Z (({ξ ξ}_{11} - - 0.5 0.5)) \\ . . \\ . . \\ . . \\ Z Z (({ξ ξ}_{n no} - - 0.5 0.5)) \end{matrix}] + + τ τ \cdot &Center Dot; {I I}_{n no + + 11} &CircleTimes; &CircleTimes; {B B}_{22} [\begin{matrix} Z Z (({ξ ξ}_{00} - - 11)) \\ Z Z (({ξ ξ}_{11} - - 11)) \\ . . \\ . . \\ . . \\ Z Z (({ξ ξ}_{n no} - - 11)) \end{matrix}] - - - - - - ((88))

in, Denotes the Kronecker product, ie I _n+1 means an n+1-dimensional square matrix.

6. the multi-turn tool parallel turning stability determination method based on differential quadrature method according to claim 5, is characterized in that, described step 5, is specifically:

Using the differential quadrature method, the velocity term on the left side of formula (8) is represented by the displacement term; for the convenience of expression, record for f(ξ _i ) represents the scalar expression symbol of the state displacement item;

First use (n+1) points (ξ ₀ ,f(ξ ₀ )),(ξ ₁ ,f(ξ ₁ )),…,(ξ _n ,f(ξ _n )) for Lagrangian interpolation, The result is shown in formula (9):

f f ((ξ ξ)) = = {Σ Σ}_{j j = = 00}^{n no} {l l}_{j j} ((ξ ξ)) f f (({ξ ξ}_{j j})) - - - - - - ((99))

Among them, the interpolation basis function is:

{l l}_{i i} ((ξ ξ)) = = \underset{k k &NotEqual; &NotEqual; i i}{{Π Π}_{k k = = 00}^{n no}} \frac{((ξ ξ - - {ξ ξ}_{k k}))}{(({ξ ξ}_{i i} - - {ξ ξ}_{k k}))} - - - - - - ((1010))

Deriving the Lagrangian interpolation function, and then substituting the time points of the velocity item at the left end of formula (8), the result is shown in formula (11):

Among them, the expression of H is shown in formula (12):

{H h}_{i i j j} = = \frac{{dl dl}_{j j} ((ξ ξ))}{d d ξ ξ} {| |}_{ξ ξ = = {ξ ξ}_{i i}} = = \{\begin{matrix} \frac{\underset{k k &NotEqual; &NotEqual; i i,, j j}{{Π Π}_{k k = = 00}^{n no}} (({ξ ξ}_{i i} - - {ξ ξ}_{k k}))}{\underset{k k &NotEqual; &NotEqual; j j}{{Π Π}_{k k = = 00}^{n no}} (({ξ ξ}_{j j} - - {ξ ξ}_{k k}))},, & j j &NotEqual; &NotEqual; i i \\ \underset{k k &NotEqual; &NotEqual; i i}{{Σ Σ}_{k k = = 00}^{n no}} \frac{11}{{ξ ξ}_{j j} - - {ξ ξ}_{k k}},, & j j = = i i \end{matrix} - - - - - - ((1212)) . .

7. the multi-turn tool parallel turning stability determination method based on differential quadrature method according to claim 6, is characterized in that, described step 6, is specifically:

The delay displacement term Z(ξ-0.5) is represented by a Chebyshev point of the second kind on the interval [-1,0] and [0,1]; find i such that ξ _i -0.5<0 and ξ _i+1 - 0.5>0, i=0,1,...,n; take the second kind of Chebyshev points on the intervals [-1,0] and [0,1] as the interpolation points to obtain the Lagrangian interpolation function, and then Substituting the time point of the time-delay displacement term, we get:

in,

{T T}_{i i j j}^{d d} = = \{\begin{matrix} \frac{\underset{k k &NotEqual; &NotEqual; j j}{{Π Π}_{k k = = 00}^{n no}} [[(({ξ ξ}_{i i} - - 0.5 0.5)) - - (({ξ ξ}_{k k} - - 11))]]}{\underset{k k &NotEqual; &NotEqual; j j}{{Π Π}_{k k = = 00}^{n no}} [[(({ξ ξ}_{j j} - - 11)) - - (({ξ ξ}_{k k} - - 11))]]},, & (({ξ ξ}_{i i} - - 0.5 0.5)) &NotElement; &NotElement; {{{ξ ξ}_{k k} - - 11 | | k k = = 00,, 11,, ... ...,, n no}} \\ 00,, & (({ξ ξ}_{i i} - - 0.5 0.5)) &Element; &Element; {{{ξ ξ}_{k k} - - 11 | | k k = = 00,, 11,, ... ...,, j j - - 11,, j j + + 11,, ... ...,, n no}} \\ 11,, & {ξ ξ}_{i i} - - 0.5 0.5 = = {ξ ξ}_{j j} - - 11 \end{matrix} - - - - - - ((1414))

in,

{T T}_{i i j j}^{p p} = = \{\begin{matrix} \frac{\underset{k k &NotEqual; &NotEqual; j j}{{Π Π}_{k k = = 00}^{n no}} [[(({ξ ξ}_{i i} - - 0.5 0.5)) - - {ξ ξ}_{k k}]]}{\underset{k k &NotEqual; &NotEqual; j j}{{Π Π}_{k k = = 00}^{n no}} (({ξ ξ}_{i i} - - {ξ ξ}_{k k}))},, & (({ξ ξ}_{i i} - - 0.5 0.5)) &NotElement; &NotElement; {{{ξ ξ}_{k k} - - 11 | | k k = = 00,, 11,, ... ...,, n no}} \\ 00,, & (({ξ ξ}_{i i} - - 0.5 0.5)) &Element; &Element; {{{ξ ξ}_{k k} | | k k = = 00,, 11,, ... ...,, j j - - 11,, j j + + 11,, ... ...,, n no}} \\ 11,, & {ξ ξ}_{i i} - - 0.5 0.5 = = {ξ ξ}_{j j} \end{matrix} - - - - - - ((1616))

The discretization of the delay term Z(ξ-1) on the interval [-1,0] is shown in formula (17):

It can be obtained that T=I _n+1 ,

Among them, T means

8. The multi-turn tool parallel turning stability determination method based on differential quadrature method according to claim 7, characterized in that, said step 7 is specifically:

Construct the Floquet transfer matrix:

The matrices H, T ^d , T ^p , T are recorded as

Among them, T ^d represents T ^p said

Combining the above formulas, it can be concluded that the state transition relationship between two adjacent intervals [-1,0] and [0,1] is:

P P [\begin{matrix} Z Z (({ξ ξ}_{00})) \\ Z Z (({ξ ξ}_{11})) \\ . . \\ . . \\ . . \\ Z Z (({ξ ξ}_{n no})) \end{matrix}] = = Q Q [\begin{matrix} Z Z (({ξ ξ}_{00} - - 11)) \\ Z Z (({ξ ξ}_{11} - - 11)) \\ . . \\ . . \\ . . \\ Z Z (({ξ ξ}_{n no} - - 11)) \end{matrix}] - - - - - - ((1818))

in,

In formula (19) and formula (20), 0 _4×(4n+4) represents the zero matrix whose dimension is 4×(4n+4), I _4×4 represents the identity matrix whose dimension is 4×4, 0 _4×4 means zero matrix with dimension 4×4, 0 _4n×4 means zero matrix with dimension 4n×4, I _n×n means identity matrix with dimension n×n, 0 _4n×4n means The dimension is a zero matrix of 4n×4n, and 0 _4×(4n+4) represents a zero matrix whose dimension is 4×(4n+4);

The Floquet state transition matrix Φ is shown in the formula:

in, Represents the Penrose-Moor generalized inverse;

According to Floquet theory, if the modulus of all eigenvalues of Φ are less than 1, the system is stable; if the modulus of any eigenvalue of Φ is greater than 1, the system is unstable.

9. the method for judging stability of multi-turning tool parallel turning based on differential quadrature method according to claim 1, is characterized in that, also comprises the steps:

Step 8: And draw the stability map of the system in the delay parameter space.