CN116484533A - 基于伯努利分布与混合驱动法的铣削稳定性分析方法 - Google Patents
基于伯努利分布与混合驱动法的铣削稳定性分析方法 Download PDFInfo
- Publication number
- CN116484533A CN116484533A CN202310443538.9A CN202310443538A CN116484533A CN 116484533 A CN116484533 A CN 116484533A CN 202310443538 A CN202310443538 A CN 202310443538A CN 116484533 A CN116484533 A CN 116484533A
- Authority
- CN
- China
- Prior art keywords
- milling
- stability
- state
- correction function
- sld
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000003801 milling Methods 0.000 title claims abstract description 91
- 238000000034 method Methods 0.000 title claims abstract description 54
- 238000009826 distribution Methods 0.000 title claims abstract description 30
- 238000004458 analytical method Methods 0.000 title claims abstract description 14
- 238000012937 correction Methods 0.000 claims abstract description 47
- 239000011159 matrix material Substances 0.000 claims abstract description 42
- 238000001228 spectrum Methods 0.000 claims abstract description 25
- 238000005457 optimization Methods 0.000 claims abstract description 21
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 20
- 230000002068 genetic effect Effects 0.000 claims abstract description 19
- 230000003595 spectral effect Effects 0.000 claims abstract description 19
- 230000007704 transition Effects 0.000 claims abstract description 19
- 238000005520 cutting process Methods 0.000 claims description 31
- 238000013016 damping Methods 0.000 claims description 15
- 230000000694 effects Effects 0.000 claims description 14
- 238000002474 experimental method Methods 0.000 claims description 13
- 238000013528 artificial neural network Methods 0.000 claims description 11
- 238000004088 simulation Methods 0.000 claims description 10
- 238000004364 calculation method Methods 0.000 claims description 8
- 230000008569 process Effects 0.000 claims description 6
- 238000012545 processing Methods 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 6
- 230000001687 destabilization Effects 0.000 claims description 3
- 238000006073 displacement reaction Methods 0.000 claims description 3
- 230000000737 periodic effect Effects 0.000 claims description 3
- 230000008929 regeneration Effects 0.000 claims description 3
- 238000011069 regeneration method Methods 0.000 claims description 3
- 230000006870 function Effects 0.000 description 45
- 238000010586 diagram Methods 0.000 description 13
- 210000002569 neuron Anatomy 0.000 description 7
- 230000003068 static effect Effects 0.000 description 3
- 238000012549 training Methods 0.000 description 3
- 230000004913 activation Effects 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 229910000838 Al alloy Inorganic materials 0.000 description 1
- 230000001133 acceleration Effects 0.000 description 1
- 125000003275 alpha amino acid group Chemical group 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000003754 machining Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000013433 optimization analysis Methods 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 230000001172 regenerating effect Effects 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000009827 uniform distribution Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/17—Mechanical parametric or variational design
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/04—Architecture, e.g. interconnection topology
- G06N3/048—Activation functions
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/08—Learning methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/12—Computing arrangements based on biological models using genetic models
- G06N3/126—Evolutionary algorithms, e.g. genetic algorithms or genetic programming
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/08—Probabilistic or stochastic CAD
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Evolutionary Computation (AREA)
- Computational Mathematics (AREA)
- Software Systems (AREA)
- Health & Medical Sciences (AREA)
- Biophysics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Computing Systems (AREA)
- Geometry (AREA)
- Artificial Intelligence (AREA)
- Molecular Biology (AREA)
- Computational Linguistics (AREA)
- Biomedical Technology (AREA)
- General Health & Medical Sciences (AREA)
- Bioinformatics & Computational Biology (AREA)
- Computer Hardware Design (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Evolutionary Biology (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Medical Informatics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physiology (AREA)
- Genetics & Genomics (AREA)
- Operations Research (AREA)
- Feedback Control In General (AREA)
Abstract
本发明公开了一种基于伯努利分布与混合驱动法的铣削稳定性分析方法，建立铣削动力学模型构造铣削稳定性的判断模型，建立伯努利分布的概率函数判断铣削过程的稳定性，构造修正函数，建立谱半径的数据驱动模型，遗传算法优化其它属性参数，最终得到优化后的属性参数进而绘制SLD图。本发明通过理论分析使用状态转移矩阵的谱半径进行铣削稳定性判断。通过分析SLD图的性质，使用两种不同的优化方法先快速确定固有圆频率绘制准确性较高的SLD图，再综合考虑准确性和稳定性确定其它属性参数。
Description
技术领域
本发明涉及铣削稳定性叶瓣图优化分析技术，特别涉及一种基于伯努利分布与混合驱动法的铣削稳定性分析方法。
背景技术
由于大型航空整体结构件尺寸大、材料切除率高、薄壁部位多、刚性差，在高速铣削过程中不仅容易引起加工变形，而且极易发生振动，导致加工精度和加工表面质量难以控制。由再生效应引起的切削厚度的动态变化是铣削颤振的最重要因素，对于铣削加工过程中不可避免产生的颤振问题，较为广泛应用以指导实际加工的方法是通过稳定性叶瓣图进行铣削稳定性分析。
目前，相关研究主要探讨如何提高SLD图的准确性。由于基于物理模型方法存在输入参数不准确的现象，而基于数据驱动方法缺乏通用性和物理可解释性，鲜有综合两种方法的特点进行铣削稳定性的分析。
发明内容
为了准确进行铣削稳定性分析，找到合理的颤振和稳定的边界条件，本发明结合物理模型方法和数据驱动方法的优势提出一种混合驱动的分析方法进行铣削稳定性的预测，最终进行稳定性叶瓣图的重绘找到颤振稳和定边界。
本发明采用以下技术方案实现上述目的。基于伯努利分布与混合驱动法的铣削稳定性分析方法，其步骤如下：
步骤一：建立铣削动力学模型构造铣削稳定性的判断模型，进而使用谱半径进行颤振稳定性的求解，
S01：通过二阶时滞微分方程描述考虑再生效应的动态铣削过程：
式中，M是模态质量矩阵，q(t)是刀具振动位移矢量，C是模态阻尼矩阵，K是模态刚度矩阵，ap是轴向切削深度，D(t)是动态铣削力矩阵，T是一个刀齿的切削周期；
S02：使用状态空间形式方程描述动态铣削过程：
式中，x(t)＝[q(t),p(t)]T，A是仅与模态参数有关的常系数矩阵，B(t)是仅与动态铣削力有关的周期系数矩阵；
S03：描述任一时间区间内的振动方程：
式中，ti是离散时间，且不切削时的振动方程为
S04：分别表示t1,t2,t3,t4时刻的状态项x1,x2,x3,x4：
式中，tf是初始时刻到刀齿开始切削的时间段；
S05：进行状态项和时滞项的分离：
离散点ti+4处的状态项xi+4表示：
分离式(5)的状态项和时滞项：
式中：
S06：连续形式状态方程的离散形式：
联立式(4,6)，表示式(2)中连续形式状态方程的离散形式：
GXm+1＝HXm+1-T (7)
一个铣削周期内状态转移矩阵的表示：
Ψ＝G-1H (8)
式(8)是铣削稳定性判断模型；稳定性根据状态转移矩阵Ψ的谱半径进行判断；
表1状态转移矩阵谱半径与系统状态的关系
| 状态转移矩阵谱半径 | 系统状态 |
| λ(Ψ)＜1 | 系统处于稳定状态 |
| λ(Ψ)＞1 | 系统处于失稳状态 |
| λ(Ψ)＝1 | 系统处于临界状态 |
步骤二：建立伯努利分布的概率函数判断铣削过程的稳定性：由于实际情况下实验参数采集不完全准确，则通过步骤一计算谱半径判断铣削状态不准确，因此通过实验结合稳定性判断模型修正实验参数提高预测准确性，具体过程如下：
S01：通过谱半径进行铣削稳定性结果的判断：
通过实验求解状态转移矩阵Ψ的谱半径需要铣削工艺参数x＝(Ω,ap)和属性参数θ＝(ωx,ωy,cx,cy,kx,ky,Dt,Dr)，其中：Ω是主轴转速，ap是铣削深度，ωx,ωy是固有圆频率，cx,cy是阻尼比，kx,ky是模态刚度，Dt,Dr是切削力系数；
假设根据式(8)可知，铣削稳定性结果如下：
式中，y是稳定性状态，λ(Ψ)是状态转移矩阵Ψ的谱半径；
S02：将稳定性结果离散量等效转化为连续量，通过概率形式进行铣削稳定性表示：
式中，p(y|θ,x)是铣削稳定性判断结果，是颤振概率，颤振概率的数值是基于谱半径性质构造的逻辑函数如式(11)：
式中，η是超参数；
S03：通过伯努利分布规律进行预测结果的判断：
λ(θ,x)＝1时将0.5作为临界值进行如表2的判断：无论实际情况如何，若稳定性判断正确，则p(y|θ,x)的数值恒大于0.5；
表2伯努利分布规律
步骤三：构造修正函数：进行n次后续颤振实验，更新稳定性叶瓣图，获得更为合理的属性参数，具体步骤如下：
S01：构造预修正函数f(θ,x)：
进行n次颤振实验，通过伯努利分布构造预修正函数：
式中，u(ζ)是判断函数，f(θ,x)是预修正函数；
若预测结果与实验结果相同则判断正确f(θ,x)增加1，即预修正函数数值与判断正确的实验组数值相同；
S02：构造再修正函数g(θ,x)：
预修正函数的不足之处，在预修正曲线1和2之间，可能存在不同铣削工况对应的颤振点，为了获得安全性更高的SLD，从加工精度和效率的角度，需要保证在SLD曲线下方尽可能稳定；
铣削条件分布点数值越大，SLD图越稳定，构造再修正函数：
式中，v(ζ)是判断函数，g(θ,x)是再修正函数；
S03：构造安全性修正函数h(θ,x)：
f(θ,x)表示判断正确的个数，g(θ,x)表示判断正确的安全性，构造安全性修正函数：
h(θ,x)＝f(θ,x)+g(θ,x) (16)
式中，h(θ,x)是安全性修正函数；
步骤四：建立谱半径的数据驱动模型：提高计算速度并推广到所有铣削工艺参数的区间范围，具体步骤如下：
S01：在误差范围内随机采样属性参数进行模拟实验，建立与谱半径之间的联系；
S02：确定BP神经网络拓扑结构；
S03：建立BP神经网络的拓扑结构后，设置相关参数进行网络模型的训练；
步骤五：黄金分割法优化固有圆频率：
误差范围令f(θ,x)最大，采用黄金分割法对固有圆频率进行优化调整至理想状态，
式中，σω是固有圆频率允许误差；
S01：设置黄金分割法的上下边界ωL、ωR；
S02：设置每次通过黄金分割法更新后可能作为新误差范围的上、下边界ω′L、ω′R；
S03：设置黄金分割法的更新方向：
每次更新可能的方向函数direct(·)定义如下：在[ωL,(ωL+ωR)/2]进行10次等距采样，每次采样计算一个f(θi,x)得到同样在[(ωL+ωR)/2,ωR]进行10次等距采样得到/>选择数值较大的一方作为更新方向；
S04：设置阈值进行更新优化，(ωL+ωR)/2的结果作优化后的最佳固有圆频率；
步骤6：遗传算法优化其它属性参数：
误差范围令h(θ,x)最大，采用遗传算法对阻尼比、模态刚度和切削力系数进行优化调整至理想状态；
式中，π是其它属性参数，σπ是其它属性参数允许误差；
S01：设置遗传算法的参数：
S02：使用遗传算法进行优化，适应度定义为：
式中，e(U)是适应度；
最终得到优化后的属性参数θ进而绘制SLD图。
本发明通过理论分析使用状态转移矩阵的谱半径进行铣削稳定性判断。通过概率形式将稳定性分析的离散变量转化为连续的概率变量便于优化，且通过伯努利分布规律对于实际情况做出正确的判断结果。在构造修正函数的同时，不仅考虑预测的准确性，且保证了预测的安全性。通过建立谱半径的数据驱动模型，使用误差构造模拟实验扩大预测模型的输入样本，相较于一般的物理模型提高计算速度，并推广到所有属性参数的连续区间范围，即在确定铣削工艺参数且进行试验后能够快速求解该条件下对应的SLD图进行稳定性判断。通过分析SLD图像的性质，使用两种不同的优化方法先快速确定固有圆频率绘制准确性较高的SLD图，再综合考虑准确性和稳定性确定其它属性参数。
附图说明
图1是本发明的铣削动力学模型图；图2是本发明的颤振概率函数图；
图3是本发明中相同的预修正函数值示意图；
图4a是颤振概率大于0.5的SLD曲线；图4b是颤振概率大于0.5的SLD曲线；
图5是网络拓扑结构图；
图6a是固有圆频率对SLD的影响；图6b是阻尼比对SLD的影响；
图6c是模态刚度对SLD的影响；图6d是切削力系数对SLD的影响；
图7是进行动态优化的数据点；
图8a是优化前是否颤振的SLD图；图8b是优化前的SLD图预测准确率；
图9a是黄金分割法优化后是否颤振的SLD图；
图9b是黄金分割法优化后的SLD图预测准确率；
图10a是遗传算法优化后是否颤振的SLD图；
图10b是遗传算法优化后的SLD图预测准确率；
具体实施方式
以下结合附图和实施例对本发明做进一步说明。
一种基于伯努利分布与混合驱动法的铣削稳定性分析方法，其步骤如下：
步骤一：建立铣削动力学模型(如图1所示)，F是受到的铣削力，fz是进给速度，Ω是主轴转速，k是模态刚度，c是阻尼比。
S01：考虑再生效应的动态铣削过程：
通过如式(1)的二阶时滞微分方程来描述：
式(1)中，M是模态质量矩阵，C是模态阻尼矩阵，K是模态刚度矩阵，q(t)是刀具振动位移矢量，f(t)是静态铣削力矩阵，D(t)是动态铣削力矩阵，ap是轴向切削深度，T＝60/NΩ是一个刀齿的切削周期，N为刀齿个数。
由于静态力不影响系统的稳定性，式(1)可退化为省略静态铣削力矩阵的形式：
式(2)中，模态质量矩阵M、模态刚度矩阵K与固有圆频率ω之间的关系为：
|K-ω2M|＝0 (3)
S02：将二阶时滞微分方程转换为状态空间形式方程：
若令x(t)＝[q(t),p(t)]T，则式(2)可表示为：
式(4)中，是仅与模态参数有关的常系数矩阵，是仅与动态铣削力有关的周期系数矩阵，且满足条件B(t)＝B(t+T)。
S03：将连续时间表示的方程转化为离散时间表示的方程：
刀齿上一周期切削结束离开工件时刻为当前周期的初始时刻t0，从t0到刀齿开始切削的时间段为tf，在t＝t0+tf时切入工件，直到当前周期结束的时间段为T-tf。在tf时间段内刀齿自由振动，在T-tf时间段内刀齿强迫振动。一个刀齿周期T可划分为非切削时间段和切削时间段，将强迫振动时间段离散为m个区间长度为τ的时间段，并将连续时间表示为离散时间如式(5)：
ti＝t0+tf+(i-1)τ (5)
式中，i＝1,2,...,m,m+1。
将式(4)中的apB(t)[x(t)-x(t-T)]看作齐次方程的非齐次项，转化为式(6)：
若t∈[ti,ti+1]，由式(6)描述任一时间区间内振动方程如式(7)：
由于在[t0,t0+tf]内刀具不切削，故B(s)＝0，式(6)退化为式(8)：
S04：表示t1时刻的状态项x1：
记xi＝x(ti)，xi-T＝x(ti-T)，Bi＝B(ti)，Bi-T＝B(ti-T)，将式(8)表示为式(9)：
S05：表示t2时刻的状态项x2：
在区间[t1,t2]上，由式(7)将离散点t2处的状态项x2表示为式(10)：
式(10)中通过梯形求积公式进行近似为式(11)：
式(11)根据状态项和时滞项进行归类整理如式(12)：
S06：表示t3时刻的状态项x3：
类似地，由式(7)将离散点t3处的状态项x3示为式(13)：
式(13)中通过Simpson求积公式近似为式(14)：
同样，式(14)根据状态项和时滞项进行归类整理如式(15)：
S07：表示t4时刻的状态项x4：
类似地，由式(7)将离散点t4处的状态项x4示为式(16)：
式(16)通过Newton求积公式进一步描述如式(17)：
在区间[ti,ti+4]上，由式(7)将离散点ti+4处的状态项xi+4表示为式(18)：
S08：由Cotes求积公式分离状态项和时滞项：
分离式(18)的状态项和时滞项如下：
式(19)中：
S09：求解最终连续形式的状态方程的离散形式：
联立式(9,12,15,17,19)，计算式(4)连续形式状态方程的离散形式为式(22)：
GXm+1＝HXm+1-T (22)
式(22)中
因此，一个铣削周期内的状态转移矩阵可以表示为：
Ψ＝G-1H (26)
由Floquet理论，系统稳定性可以根据状态转移矩阵Ψ的谱半径大小来判断(如表1所示)，式(26)被称之为铣削稳定性的判断模型。
表1状态转移矩阵谱半径与系统状态的关系
| 状态转移矩阵谱半径 | 系统状态 |
| λ(Ψ)＜1 | 系统处于稳定状态 |
| λ(Ψ)＞1 | 系统处于失稳状态 |
| λ(Ψ)＝1 | 系统处于临界状态 |
步骤二：建立基于伯努利分布的概率函数。分析铣削过程的稳定性，同时需要铣削工艺参数x＝(Ω,ap)和属性参数θ＝(ωx,ωy,cx,cy,kx,ky,Dt,Dr)。由于铣削过程中受多种因素相互制约而存在测量数据不准确问题，需要对属性参数进行修正。假设属性参数误差σ＝(σωx,σωy,σcx,σcy,σkx,σky,σDt,σDr)。建立伯努利分布的概率函数能够应用于误差范围内属性参数的修正，具体过程如下：
S01：通过谱半径进行铣削稳定性结果的判断：
通过式(26)将铣削稳定性结果进行式(27)的表示：
式(27)中，y是稳定性状态，其中0表示稳定，1表示失稳，λ(Ψ)是状态转移矩阵Ψ的谱半径。
S02：将铣削稳定性结果通过概率形式进行表示：
连续量优化在动态优化问题中具有更好的效果。利用伯努利分布将稳定性结果y中0或1的离散量等效转化为概率值的连续量。定义/>为颤振的概率，/>为稳定的概率，将式(27)中铣削稳定性结果用概率形式表示为式(28)：
S03：构造谱半径与颤振概率的函数关系：
颤振概率合适的表示方法是实现稳定性结果y与概率p之间等效表示的关键因素。基于谱半径的性质构造逻辑函数如式(29)，超参数η的取值会影响目标趋近的速度(如附图2所示)，这里取η＝4：
S04：通过伯努利分布规律进行预测结果的判断：
λ(θ,x)＝1时将0.5作为临界值进行如下判断：无论实际实验结果如何，只要稳定性判断正确，p(y|θ,x)的数值恒大于0.5(如表2所示)。表2伯努利分布规律
步骤三：构造修正函数，进行n次后续颤振实验，用于更新修正稳定性叶瓣图，进而获得更为合理的属性参数。具体步骤如下：
S01：构造预修正函数f(θ,x)：
设置阶跃函数u(ζ)作为p(yi|θ,xi)的激活函数，保留正确的判断结果并输出1，过滤错误的判断结果。进行n次颤振实验获得伯努利分布数值，将其通过激活函数u(ζ)处理再求和构造预修正函数f(θ,x)，如式(30,31)：
若用于动态修正的实验有n组，稳定性判断正确的实验有m组，则模型判断的准确率为m/n。每完成一次正确的判断，f(θ,x)增加1，预修正函数数值与判断正确的实验组数值相同，模型判断的准确率也可以表示为f(θ,x)/n。
S02：构造再修正函数g(θ,x)：
预修正函数具有如下不足之处(如图3所示)，图中包含两条预修正的SLD曲线1和曲线2，圆圈表示稳定，方框表示颤振，两条曲线均作出正确判断，所对应的预修正函数与/>的值相等，无法判断其修正效果的好坏。在预修正曲线1和2之间，可能存在不同铣削工况对应的颤振点，为了获得安全性更高的SLD，从加工精度和效率的角度，需要保证在SLD曲线下方尽可能保持稳定。选择图中的预修正曲线1作为最终的修正SLD曲线，效果要强于预修正曲线2。
每个铣削条件分布点所对应的稳定性状态已转化0～1之间颤振概率的数值，概率＞0.5的图像(如图4a所示)，概率＜0.5的图像(如图4b所示)。对于铣削条件分布点，/>的数值越大，SLD图的结果越稳定。进一步构造出再修正函数g(θ,x)如下：
/>
式中v(ζ)是ReLU函数，可仅保留判断正确部分的具体数值。若颤振即yi＝1，取较大值，v(ζ)p(yi|θ,xi)也取较大值，则g(θ,x)取较大值；若颤振即yi＝0，/>取较大值，v(ζ)p(yi|θ,xi)取较小值，1-v(ζ)p(yi|θ,xi)取较大值，则g(θ,x)仍取较大值。
S03：构造安全性修正函数h(θ,x)：
从修正函数构造的过程来看，f(θ,x)表示判断正确的个数，g(θ,x)表示判断正确的安全性，且g(θ,x)的数值不影响判断正确的个数两者彼此相互不影响。可以构造出更为合理的SLD安全性修正函数h(θ,x)如下：
h(θ,x)＝f(θ,x)+g(θ,x) (34)
步骤四：在更新计算修正函数的过程中，每次迭代都需要重新计算谱半径λ。若每次更新均采用式(26)将花费较长的计算时间，为此建立谱半径的数据驱动模型，计算速度远优于基于物理模型法的计算速度。具体步骤如下：
S01：通过误差进行模拟实验：
对于给定的铣削工况，在误差范围内，随机采样属性参数，利用式(26)进行谱半径模拟实验后，选择BP神经网络建立属性参数与谱半径之间的联系，即建立一定铣削工况下的谱半径神经网络预测模型。
S02：确定BP神经网络拓扑结构：
网络结构包括1个输入层、n个隐藏层和1个输出层，输入层有8个神经元θi(1≤i≤8)，输出层只有1个神经元，即需要预测的谱半径Λ(如图5所示)。
第一层隐含层的神经元数目m1按照通用的经验公式确定为m1＝17。从第二层隐藏层开始，神经元的个数通过网络结构增长方法来确定，即依次增加神经元的个数，直到预测结果误差和实际误差无明显增大后确定最终的网络拓扑结构。
S03：进行网络模型的训练：
建立BP神经网络的拓扑结构后，设置相关参数进行网络模型的训练。
步骤五：黄金分割法优化固有圆频率：
误差范围内，调整固有圆频率参数至理想状态，使得f(θ,x)最大从而修正SLD。由于SLD图的左右位置仅与固有圆频率有关(如图6a所示)，故对于式(35)，采用黄金分割法对固有圆频率进行优化。
S01：设置黄金分割法的上下边界：
设置固有圆频率ω在初始误差σω范围的上、下边界ωL、ωR。
S02：设置黄金分割法的更新边界：
设置每次通过黄金分割法更新后可能作为新误差范围的上、下边界ω′L、ω′R。
S03：设置黄金分割法的更新方向：
每次更新可能的方向函数direct(·)定义如下：在[ωL,(ωL+ωR)/2]区间上进行10次等距采样，每次采样计算一个f(θi,x)，得到同样在[(ωL+ωR)/2,ωR]区间上也进行10次等距采样，得到/>选择数值较大的作为更新方向。
S04：设置阈值进行更新优化：
迭代过程直至满足当前圆频率误差不超过阈值εω这一终止条件。(ωL+ωR)/2的结果作为黄金分割法优化后的固有圆频率的最佳数值。
步骤6：遗传算法优化其它属性参数：
误差范围内，调整阻尼比、模态刚度和切削力系数参数至理想状态，使得h(θ,x)最大从而修正SLD。阻尼比、模态刚度和切削力系数对SLD图的影响(如图6b、图6c和图6d所示)，对于式(36)，采用遗传算法对阻尼比、模态刚度和切削力系数进行优化。
S01：设置遗传算法的参数：
设置种群规模，变量的取值上下限，适应度函数，选择函数，交叉概率和变异概率。
S02：使用遗传算法进行优化：
采用遗传算法求解式(37)，找到使目标函数h(θ,x)最大时所对应的阻尼比、模态刚度和切削力系数。适应度可定义为：
实施例：选择由铣床NCT EmR-610Ms，螺旋角β＝30°、直径为16mm的两齿立铣刀，以及铝合金2024-T351工件组成的铣削系统进行铣削实验，为便于观察效果仅优化SLD图的单个叶瓣图，验证所提出的基于伯努利分布和混合驱动法的铣削稳定性分析方法。
步骤1：属性参数分析：
切削力由Kistler 9129AA多分量测力仪、NI-9234采集卡采集，模态参数由Endevco力锤2302-10、PCB加速度传感器352C32获得，实验测得的属性参数，见表3。
表3属性参数的实验值
S01：利用属性参数通过Cotes积分法获得SLD图：
利用如表3所示的铣削系统属性参数，通过式(26)Cotes积分法获得SLD图，见图6a、图6b、图6c和图6d中的黑色实线。
S02：分析固有圆频率对SLD图的影响：
固有圆频率ωx,ωy对SLD图的影响(如图6a所示)，图例中90％ωx＝4426rad/s，90％ωy＝4257rad/s，依次类推。由图可知固有圆频率主要改变SLD图的左右位置。
S03：分析除固有圆频率外其它因素对SLD图的影响：
阻尼比、模态刚度和切削力系数对SLD图的影响(如图6b、图6c、图6d所示)，由图可知除固有圆频率外其它因素的数值变化主要改变SLD图的高度和宽度。
步骤2：谱半径计算：
S01：通过属性参数及其误差进行模拟实验：
在误差范围内，对属性参数采用均匀分布的随机数方法进行1000次模拟实验。
S02：求解各组属性参数条件下不同铣削工况所对应的谱半径：
图7中是进行动态优化数据点的散点，圆圈代表稳定，叉号代表颤振，三角代表无法确定结果。设置主轴转速的网格节点步长为100，切削深度的网格节点步长为0.125，横坐标为主轴转速离散为85个节点，纵坐标为切削深度离散为27个节点，坐标平面中包含2295个网格节点，每个网格节点对应一种铣削工况，通过1000次模拟实验求解1000组属性参数条件下2295种铣削工况所对应的谱半径。
S03：通过BP神经网络建立谱半径的预测模型拓扑结构：
1000次模拟实验中每次属性参数的SLD图包含2295个网格节点，建立1000*2295个神经网络预测模型。输入8个神经元即属性参数，输出1个神经元即谱半径，设置网络结构为8-17-8-1。
S04：设置BP神经网络相关参数进行模型的训练：
设置网络激活函数分别为F1＝logsig、F2＝logsig和F3＝purelin，最大迭代周期为100次，学习率为0.01；对输入参数进行标准化预处理，选择训练算法trainlm。随机从1000组抽样模拟实验获得的950组作为训练集，剩下的50组作为测试集。
步骤3：黄金分割法动态优化单叶瓣图：
见图7,这里仅使用稳定和颤振的数据点进行更新和修正SLD图。选择铣削工况的网络分布恰好覆盖SLD的一个叶瓣的43次铣削加工实验，优化前的SLD图(如图8a所示)，圆圈代表稳定，叉号代表颤振，其准确率见图8b所示，下三角代表预测正确，上三角代表预测错误，其它附图依次类推，化前的SLD准确性是70％。黄金分割法优化后的SLD图(如图9a所示)，设定阈值εω为1e-6，其准确性是93％(如图9b所示)。
步骤4：遗传算法动态优化单叶瓣图：
保持固有圆频率的最优值不变，进一步对其它6个属性参数进行优化。遗传算法优化后的SLD图(如图10a所示)，其准确率(如图10b所示)，设置种群规模为50，变量的取值上下限设置为±20％，输入参数为两个方向的阻尼比、模态刚度和切削力系数参数，优化目标为使得适应度函数的数值即修正函数的数值最大，选择函数为随机均匀选择，交叉概率为0.8，变异概率为0.2。通过遗传算法优化后的SLD准确性是98％。所有模型参数经过2次优化后的数值见表4。
表4优化单瓣后的属性参数
/>
Claims (1)
1.基于伯努利分布与混合驱动法的铣削稳定性分析方法，其特征在于，步骤如下：
步骤一：建立铣削动力学模型构造铣削稳定性的判断模型，进而使用谱半径进行颤振稳定性的求解，
S01：通过二阶时滞微分方程描述考虑再生效应的动态铣削过程：
式中，M是模态质量矩阵，q(t)是刀具振动位移矢量，C是模态阻尼矩阵，K是模态刚度矩阵，ap是轴向切削深度，D(t)是动态铣削力矩阵，T是一个刀齿的切削周期；
S02：使用状态空间形式方程描述动态铣削过程：
式中，x(t)＝[q(t),p(t)]T，A是仅与模态参数有关的常系数矩阵，B(t)是仅与动态铣削力有关的周期系数矩阵；
S03：描述任一时间区间内的振动方程：
式中，ti是离散时间，且不切削时的振动方程为
S04：分别表示t1,t2,t3,t4时刻的状态项x1,x2,x3,x4：
式中，tf是初始时刻到刀齿开始切削的时间段；
S05：进行状态项和时滞项的分离：
离散点ti+4处的状态项xi+4表示：
分离式(5)的状态项和时滞项：
式中：
S06：连续形式状态方程的离散形式：
联立式(4,6)，表示式(2)中连续形式状态方程的离散形式：
GXm+1＝HXm+1-T (7)
一个铣削周期内状态转移矩阵的表示：
Ψ＝G-1H (8)
式(8)是铣削稳定性判断模型；稳定性根据状态转移矩阵Ψ的谱半径进行判断；
表1状态转移矩阵谱半径与系统状态的关系
步骤二：建立伯努利分布的概率函数判断铣削过程的稳定性：由于实际情况下实验参数采集不完全准确，则通过步骤一计算谱半径判断铣削状态不准确，因此通过实验结合稳定性判断模型修正实验参数提高预测准确性，具体过程如下：
S01：通过谱半径进行铣削稳定性结果的判断：
通过实验求解状态转移矩阵Ψ的谱半径需要铣削工艺参数x＝(Ω,ap)和属性参数θ＝(ωx,ωy,cx,cy,kx,ky,Dt,Dr)，其中：Ω是主轴转速，ap是铣削深度，ωx,ωy是固有圆频率，cx,cy是阻尼比，kx,ky是模态刚度，Dt,Dr是切削力系数；
假设根据式(8)可知，铣削稳定性结果如下：
式中，y是稳定性状态，λ(Ψ)是状态转移矩阵Ψ的谱半径；
S02：将稳定性结果离散量等效转化为连续量，通过概率形式进行铣削稳定性表示：
式中，p(y|θ,x)是铣削稳定性判断结果，是颤振概率，颤振概率的数值是基于谱半径性质构造的逻辑函数如式(11)：
式中，η是超参数；
S03：通过伯努利分布规律进行预测结果的判断：
λ(θ,x)＝1时将0.5作为临界值进行如表2的判断：无论实际情况如何，若稳定性判断正确，则p(y|θ,x)的数值恒大于0.5；
表2伯努利分布规律
步骤三：构造修正函数；进行n次后续颤振实验，更新稳定性叶瓣图，获得更为合理的属性参数，具体步骤如下：
S01：构造预修正函数f(θ,x)：
进行n次颤振实验，通过伯努利分布构造预修正函数：
式中，u(ζ)是判断函数，f(θ,x)是预修正函数；
若预测结果与实验结果相同则判断正确f(θ,x)增加1，即预修正函数数值与判断正确的实验组数值相同；
S02：构造再修正函数g(θ,x)：
预修正函数的不足之处，在预修正曲线1和2之间，可能存在不同铣削工况对应的颤振点，为了获得安全性更高的SLD，从加工精度和效率的角度，需要保证在SLD曲线下方尽可能稳定；
铣削条件分布点数值越大，SLD图越稳定，构造再修正函数：
式中，v(ζ)是判断函数，g(θ,x)是再修正函数；
S03：构造安全性修正函数h(θ,x)：
f(θ,x)表示判断正确的个数，g(θ,x)表示判断正确的安全性，构造安全性修正函数：
h(θ,x)＝f(θ,x)+g(θ,x) (16)
式中，h(θ,x)是安全性修正函数；
步骤四：建立谱半径的数据驱动模型：提高计算速度并推广到所有铣削工艺参数的区间范围，具体步骤如下：
S01：在误差范围内随机采样属性参数进行模拟实验，建立与谱半径之间的联系；
S02：确定BP神经网络拓扑结构；
S03：建立BP神经网络的拓扑结构后，设置相关参数进行网络模型的训练；
步骤五：黄金分割法优化固有圆频率：
误差范围令f(θ,x)最大，采用黄金分割法对固有圆频率进行优化调整至理想状态，
式中，σω是固有圆频率允许误差；
S01：设置黄金分割法的上下边界ωL、ωR；
S02：设置每次通过黄金分割法更新后可能作为新误差范围的上、下边界ω′L、ω′R；
S03：设置黄金分割法的更新方向：
每次更新可能的方向函数direct(·)定义如下：在[ωL,(ωL+ωR)/2]进行10次等距采样，每次采样计算一个f(θi,x)得到同样在[(ωL+ωR)/2,ωR]进行10次等距采样得到/>选择数值较大的一方作为更新方向；
S04：设置阈值进行更新优化，(ωL+ωR)/2的结果作优化后的最佳固有圆频率；
步骤6：遗传算法优化其它属性参数：
误差范围令h(θ,x)最大，采用遗传算法对阻尼比、模态刚度和切削力系数进行优化调整至理想状态；
式中，π是其它属性参数，σπ是其它属性参数允许误差；
S01：设置遗传算法的参数：
S02：使用遗传算法进行优化，适应度定义为：
式中，e(U)是适应度；
最终得到优化后的属性参数θ进而绘制SLD图。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310443538.9A CN116484533A (zh) | 2023-04-24 | 2023-04-24 | 基于伯努利分布与混合驱动法的铣削稳定性分析方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310443538.9A CN116484533A (zh) | 2023-04-24 | 2023-04-24 | 基于伯努利分布与混合驱动法的铣削稳定性分析方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116484533A true CN116484533A (zh) | 2023-07-25 |
Family
ID=87211440
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310443538.9A Pending CN116484533A (zh) | 2023-04-24 | 2023-04-24 | 基于伯努利分布与混合驱动法的铣削稳定性分析方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116484533A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117270457A (zh) * | 2023-09-26 | 2023-12-22 | 北京航空航天大学 | 数据机理混合驱动的机器人铣削稳定性建模方法 |
-
2023
- 2023-04-24 CN CN202310443538.9A patent/CN116484533A/zh active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117270457A (zh) * | 2023-09-26 | 2023-12-22 | 北京航空航天大学 | 数据机理混合驱动的机器人铣削稳定性建模方法 |
| CN117270457B (zh) * | 2023-09-26 | 2024-04-19 | 北京航空航天大学 | 数据机理混合驱动的机器人铣削稳定性建模方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN116484533A (zh) | 基于伯努利分布与混合驱动法的铣削稳定性分析方法 | |
| CN109977464B (zh) | 一种基于bp神经网络的活塞切削加工变形量的预测方法 | |
| CN114384861A (zh) | 一种基于多层感知机的数控机床多工况切削参数优化方法 | |
| CN105488282A (zh) | 一种基于动态加工特征的切削参数分段和变切深优化方法 | |
| US20210312262A1 (en) | Method for predicting status of machining operation | |
| WO2021174518A1 (zh) | 一种无颤振铣削加工表面形貌仿真方法 | |
| CN110221580B (zh) | 一种基于主轴数据仿真的进给速度优化方法 | |
| CN107248047B (zh) | 基于加工过程状态熵动态计算的加工过程状态评价方法 | |
| CN105446264B (zh) | 基于特征的机床精度优化设计方法 | |
| CN114492198A (zh) | 一种基于改进pso算法辅助svm算法的切削力预测方法 | |
| CN113239461B (zh) | 一种非对称结构复杂薄壁零件变形控制方法 | |
| CN111723440A (zh) | 一种薄壁件加工精度预测混合建模方法 | |
| CN113369994A (zh) | 一种高速铣削过程刀具状态监测方法 | |
| CN111975453B (zh) | 一种数值仿真驱动的加工过程刀具状态监测方法 | |
| Arapoğlu et al. | An ANN-based method to predict surface roughness in turning operations | |
| CN116560301A (zh) | 一种基于梯度优化的机床进给系统数理模型参数辨识方法 | |
| Zhan et al. | Optimal pitch angles determination of ball-end cutter for improving five-axis milling stability | |
| CN114509991A (zh) | 考虑参数不确定的数控机床切削稳定性预测与优化方法 | |
| CN113987698A (zh) | 一种基于数据驱动的机床功耗模型建模方法 | |
| CN111948977B (zh) | 一种用于不锈钢加工的多目标优化方法和系统 | |
| CN109940462B (zh) | 铣刀切削振动变化特性的检测与高斯过程模型构建方法 | |
| CN108520117A (zh) | 一种利用全离散法获取稳定性叶瓣图的方法 | |
| CN109933940B (zh) | 基于滚刀主轴振动响应模型的滚齿工艺参数优化方法 | |
| CN109446721B (zh) | 基于标识符软件线程执行顺序排列的机床工艺交互算法 | |
| CN114626015B (zh) | 一种基于高斯过程回归的薄壁结构切削颤振预测方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |