CN114626015A - 一种基于高斯过程回归的薄壁结构切削颤振预测方法 - Google Patents

一种基于高斯过程回归的薄壁结构切削颤振预测方法 Download PDF

Info

Publication number
CN114626015A
CN114626015A CN202210126159.2A CN202210126159A CN114626015A CN 114626015 A CN114626015 A CN 114626015A CN 202210126159 A CN202210126159 A CN 202210126159A CN 114626015 A CN114626015 A CN 114626015A
Authority
CN
China
Prior art keywords
model
gaussian process
ipw
thin
gpr
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.)
Granted
Application number
CN202210126159.2A
Other languages
English (en)
Other versions
CN114626015B (zh
Inventor
杨昀
杨阳
刘华辰
万敏
张卫红
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202210126159.2A priority Critical patent/CN114626015B/zh
Publication of CN114626015A publication Critical patent/CN114626015A/zh
Application granted granted Critical
Publication of CN114626015B publication Critical patent/CN114626015B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • 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
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Optimization (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Physics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Algebra (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Numerical Control (AREA)

Abstract

本发明提出了一种基于高斯过程回归的薄壁结构切削颤振预测方法，通过将高斯过程回归(GPR)与适当的正交分解相结合，提出了薄壁结构IPW动力学参数的替代模型。GPR方法用于学习一组已知的IPW动力学参数与相应工具位置之间的映射。在建立IPW模态振型的GPR模型之前，采用适当的正交分解来降低不同工具位置的模态振型向量组装的矩阵的阶数。所提出方法的计算时间主要由预测一组已知IPW动力学参数的时间和训练GPR模型的时间组成。仿真表明，本发明比现有技术消耗更少的计算时间；预测的稳定瓣图与实验结果的比较表明，本发明预测的IPW动力学参数足以预测铣削薄壁结构的稳定性。

Description

一种基于高斯过程回归的薄壁结构切削颤振预测方法
技术领域
本发明属于薄壁结构切削颤振预测领域，具体涉及一种基于高斯过程回归的薄壁结构切削颤振预测方法。
背景技术
先进飞机大量采用薄壁结构，以提高飞机的性能同时减轻重量。它们大多由多个薄壁和腔槽组成，壁高与厚度的比例和长度与厚度的比例通常分别为10:1～50:1和50:1～250:1。由于铣削工艺对于具有复杂几何形状的结构来说高效且经济，薄壁结构通常通过铣削工艺制造成最终形状。在半精铣和精铣过程中，由于这些结构的刚度低，容易产生颤振，会导致表面质量差和加工精度低。目前可采用多种不同的方法避免铣削颤振，例如在主轴或工件上安装主动/被动阻尼器抑制颤振方法，或基于颤振预测或识别技术的工艺参数优化方法。其中，基于预测稳定性叶瓣图(SLDs)选择无颤振工艺参数是最常用的方法。这种方法在过去的几十年中引起了广泛的关注。
加工中工件动力学参数对颤振预测结果的准确性有显着影响。然而，由于材料去除和刀具切削位置变化的影响，加工中工件(IPW)的动力学参数变化很大。因此，在半精加工和精加工过程中准确有效地获得变化的IPW动力学参数是一个重要的问题。
传统获取薄壁结构动力学参数主要通过实验测量手段完成。显然，对变化的IPW动力学参数进行进行重复的离线式的实验测量费时费力、不切实际。铣削薄壁结构颤振准确预测的关键问题是准确有效地获取时变IPW动力学参数。对薄壁结构的研究大多基于结构动力修改和模型降阶技术。其计算时间将随着要预测的刀具位置的数量而增加。如果只是对铣削工艺进行数值模拟，IPW动力学参数仅在刀具轨迹生成阶段确定的刀具位置获得，要预测的刀具位置数量以及计算时间仍可接受。然而，如果利用铣削工艺的动力学模型来优化轴向切削深度等加工参数，要预测的刀具位置的数量将大幅增加。对于这种情况，现有方法的计算效率仍然太低。
发明内容
本发明针对上述问题，提出了一种基于高斯过程回归的薄壁结构切削颤振预测方法。
本发明通过下述技术方案实现：一种基于高斯过程回归的薄壁结构切削颤振预测方法，该方法包括以下步骤：
首先，利用现有的方法可以获得某些刀具位置的加工中工件(IPW)动力学参数和前一刀具路径的高度e；
然后，利用特征正交分解(POD)方法对模态振型u(k)进行预处理，使用高斯过程回归方法从已知集合中学习等式F:e→ω(k)中的映射，ω(k)为固有频率；
根据材料去除和切削接触位置变化而引起的薄壁件动力学参数特性的变化，建立薄壁件铣削的动力学参数模型；所述动力学参数模型为结合高斯过程回归(GPR)和特征正交分解(POD)方法用于预测加工中工件(IPW)时变动力学参数的模型，是基于非参数核的概率模型；考虑已知的与切削位置对应的固有频率集合
Figure BDA0003500594530000021
其中ei＝[xi,zi,z′i]∈R3,
Figure BDA0003500594530000022
来自一个未知的分布；
所述动力学参数模型基于机器学习方法建立，所述机器学习方法包括：给定新的输入向量如切削位置e*，以及训练数据，线性回归模型如下:
ω(k)＝eTβ+ε
其中ε遵循一个独立的、均匀分布的高斯分布，其均值为零方差为σ2，即
Figure BDA0003500594530000023
根据数据估计了误差方差σ2和系数β；高斯过程回归(GPR)模型通过从高斯过程(GP)中引入潜在变量h(ei),i＝1,2,...,nk和显式基函数φ来解释响应；潜在变量的协方差函数获取响应的平滑度，基函数将输入的切削位置向量e投影到p维特征空间中；
高斯过程是一组随机变量，因此它们中的任何有限数量都有一个联合高斯分布；如果{h(e),e∈R3}是一个高斯过程，则给定nk个切削位置
Figure BDA0003500594530000024
随机变量
Figure BDA0003500594530000025
的联合分布是高斯分布；高斯过程由其均值函数m(e)和协方差函数κ(e,e′)定义，即如果{h(e),e∈R3}是一个高斯过程，便可得到E(h(e))＝m(e)和Cov(h(e),h(e′))＝E[{h(e)-m(e)},{h(e′)-m(e′)}]＝κ(e,e′)；
现考虑已下模型：
φ(e)Tβ+h(e),
其中
Figure BDA0003500594530000031
h(e)来自于一个协方差函数为κ(e,e′)的零均值高斯过程；φ(e)是一组基函数，它将R3中的切削位置向量e转换为Rp中的一个新的特征向量φ(e)；β表示基函数的系数为p维向量；此模型代表了一个高斯过程回归模型，响应的实例ω(k)建模为
Figure BDA0003500594530000032
高斯过程回归(GPR)模型是一个概率模型，对于每个切削位置ei都引入一个潜在变量h(ei)，使得高斯过程回归(GPR)模型非参数化；在向量形式下，这个高斯过程回归(GPR)模型等价于
Figure BDA0003500594530000033
其中
Figure BDA0003500594530000034
潜在变量
Figure BDA0003500594530000035
在高斯过程回归(GPR)模型中的联合分布如下
Figure BDA0003500594530000036
接近于线性回归模型，其中K(E,E)为
Figure BDA0003500594530000037
协方差函数κ(e,e′)也被称为核函数，它指定了两个潜在变量h(ei)和h(ei′)之间的协方差；
基于特征正交分解和高斯过程回归的振型预处理模型进行变模态振型的预测。
进一步地，所述特征正交分解(POD)方法包括：对于特定的切削位置ei，工件的第k阶模态的模态振型表示为
Figure BDA0003500594530000041
这组模态振型的平均值为
Figure BDA0003500594530000042
通过关注与平均值的偏差，形成一组新向量：
Figure BDA0003500594530000043
假设已找到一组基函数
Figure BDA0003500594530000044
系统中任何已知的模态振型都可以精确地表示为基函数的线性组合，如下所示：
Figure BDA0003500594530000045
其中
Figure BDA0003500594530000046
最后，需要预测的任何模态振型
Figure BDA0003500594530000047
都可以近似地表示为：
Figure BDA0003500594530000048
其中，
Figure BDA0003500594530000049
是保留的基函数的数量，
Figure BDA00035005945300000410
为基函数，
Figure BDA00035005945300000411
为保留的基函数个数，
Figure BDA00035005945300000412
为系数；设矩阵W(k)表示这些模态振型
Figure BDA00035005945300000413
的集合：
Figure BDA00035005945300000414
基函数
Figure BDA00035005945300000415
可以通过矩阵W(k)的奇异值分解得到：
W(k)＝L(k)Σ(k)R(k)T，
其中，
Figure BDA00035005945300000416
为左奇异矩阵，Σ(k)是奇异值矩阵，可以表示为
Figure BDA0003500594530000051
其中，
Figure BDA0003500594530000052
为第j个左奇异向量
Figure BDA0003500594530000053
的动能；
让
Figure BDA0003500594530000054
如果
Figure BDA0003500594530000055
δ≤10-3，L(k)中的第一个左奇异向量
Figure BDA0003500594530000056
获取了最大可能的动能，确定了基函数
Figure BDA0003500594530000057
和保留的基函数的数量
Figure BDA0003500594530000058
经过所述基函数
Figure BDA0003500594530000059
和保留基函数的个数
Figure BDA00035005945300000510
的确定过程，得到了一个新的映射：
Figure BDA00035005945300000511
计算出与新的切削位置e*相关的系数
Figure BDA00035005945300000512
最后，通过等式
Figure BDA00035005945300000513
计算出相应的模态振型。
进一步地，所述协方差函数κ(e,e′)由平方指数核、指数核、有理二次核或Matren类核函数定义，由一组核参数或超参数θ进行参数化。
有益的技术效果：
由于薄壁结构的动力学参数在铣削过程中瞬时变化，因此准确有效地预测加工中工件(IPW)动力学参数对于预测铣削薄壁结构的颤振稳定性至关重要。当必须预测大量切削位置的IPW动力学参数时，本发明通过将高斯过程回归(GPR)与适当的正交分解相结合，提出了薄壁结构IPW动力学参数的替代模型。GPR方法用于学习一组已知的IPW动力学参数与相应工具位置之间的映射。在建立IPW模态振型的GPR模型之前，采用适当的正交分解来降低不同工具位置的模态振型向量组装的矩阵的阶数。所提出方法的计算时间主要由预测一组已知IPW动力学参数的时间和训练GPR模型的时间组成。仿真表明，所提出的方法消耗较少的计算时间，同时，提出的方法与现有方法具有可比性。预测的稳定性叶瓣图与实验结果的比较表明，所提出的方法预测的IPW动力学参数足以预测铣削薄壁结构的稳定性。
通过将本发明的IPW动力学参数预测SLD与其他方法的预测结果进行比较，验证了本发明的准确性。与现有方法刀具动力学参数的预测和实验稳定性结果的比较如图7所示。注意，它们是图6的部分。在图7中，预测结果和测量结果在特定主轴转速下进行了比较。实验中，通过测量声音的频谱分布检测颤振。如果发生颤振，频谱应在刀具或工件的固有频率附近具有峰值，但与主轴或刀齿通过频率及其谐波不同。当声音或位移信号的频率分布在主轴或刀齿通过频率及其谐波上时，切削过程是稳定的。预测的稳定性极限与实验结果吻合良好。在图7所示的126组详细比较中，只有10组预测与测量结果不一致，一致率超过92％。这意味着所提出的方法可以准确地预测IPW动力学参数，并且可以集成到铣削薄壁工件的动力学参数模型中。
附图说明
此处所说明的附图用来提供对本发明实施例的进一步理解，构成本申请的一部分，并不构成对本发明实施例的限定。在附图中：
图1是动态模型的示意图；
图2是工件的三维立体图；
图3是使用不同自由度的有限元模型对实施例1、现有方法和实施例2中GPR模型训练程序的计算时间进行比较；
图4是本发明预测的固有频率与现有方法相关的固有频率的比较：左列为前三阶固有频率，右列为对应的相对误差；
图5是本发明预测的振型与现有方法预测的振型之间的相关性；
图6是测试工件的三维稳定性叶瓣图；
图7是主轴转速为12000r/min时沿刀具轨迹铣削曲面薄壁工件的稳定性极限。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚明白，下面结合实施例，对本发明作进一步的详细说明，本发明的示意性实施方式及其说明仅用于解释本发明，并不作为对本说明书发明的限定。
本发明提出了一种利用高斯过程回归的策略。在该策略中，高斯过程回归方法被用于学习一组已知的IPW动力学参数和相应的刀具位置之间的映射。通过现有的半解析方法或基于有限元的方法，在特定的刀具位置预测已知的IPW动力学参数。
实施例1
由于材料去除和切削接触位置变化而引起的薄壁件动力学参数的变化，在铣削过程中，刀具轨迹的位置决定了切削接触的实际位置。也就是说，切削力仅作用于图1所示坐标x和y确定的切削接触位置。此外，去除的材料由前一刀具路径z′的高度确定。在本发明中，列向量e＝[x,z,z′]T被称为切削位置，用于表示决定IPW时变动力学参数。
根据材料去除和切削接触位置变化而引起的薄壁件动力学参数的变化，建立薄壁件铣削的动力学参数模型，动力学参数模型的示意图如图1所示。由于工件和刀具刚度低，它们的偏转是不可避免的，并且它们可以在机床坐标系OXYZ中的切削力作用下在X、Y和Z方向上振动。因此，每个节点在X、Y和Z方向上都有3个自由度。工件和刀具的位移分别表示为一个3nw维向量Qw和一个3nt维向量Qt,其中，nw和nt分别表示与工件和刀具相关的节点编号。切削接触区沿刀具轴离散为na个节点。对于上一刀具路径的特定刀具位置和高度e＝[x,z,z′]T，物理空间中IPW和刀具的动力学参数方程可表示如下：
Figure BDA0003500594530000081
其中，Mc、Cc和Kc分别为部件c的3nc×3nc质量、阻尼和刚度矩阵。下标‘c’分别代表工件或刀具对应的‘w’或‘t’。‘t’表示时间，
Figure BDA0003500594530000082
和
Figure BDA0003500594530000083
分别表示节点的速度和加速度对应的3nc维向量。F表示与切削力相关的3na维向量。作用在工件和刀具上的力大小相同，但方向相反，它是使用现有方法来计算的。Dc表示实际的切削接触位置，是3nc×3na变换矩阵，其中除与实际切削接触位置相关的全一行向量外，其余行均为零行向量。
需要注意的是，每个节点的振幅和动态特性沿刀具轴线变化，在推导过程中也考虑了它们。当轴向切削深度较小时，振动幅值和动态特性沿刀具轴向的变化可以忽略不计。这样，就只能用一个节点来描述工件或刀具的振动，多接触点模型退化为单接触点模型。
在建模过程中生成了一个大型的薄壁结构。基于这种大规模动力学参数模型在物理空间中直接分析结构的振动和稳定性存在大量的计算负担，因此，物理空间中的动力学参数模型通常被转换到模态空间中。同时，由于获取商用机床的几何结构和估算主轴内部连接界面的刚度和阻尼特性的复杂性，很难精准地确定由质量、阻尼和刚度矩阵表示的刀具主轴系统的动力学参数，即式(1)中的Mt，Ct和Kt。然而，基于投影的思想，通过实验模态分析可以很容易地获得模态空间的动力学参数。最后，通过解析模态分析得到的工件模态矩阵Uw和实验模态分析得到的刀具模态矩阵Ut用于转换等式(1)进入模态空间：
Qc(t)＝UcΓc(t)(c＝t,w) (2)
其中，Γc表示mc维的模态位移向量。Uc是3nc×mc质量归一化模态矩阵，mc表示mc＜＜3nc的模式数。
通过将等式(2)带入等式(1)然后乘以Uc的转置矩阵，等式(1)转换为模态空间，如下所示：
Figure BDA0003500594530000091
其中
Figure BDA0003500594530000092
ωc(e)、ζc(e)和Ic分别为mc×mc对角固有频率矩阵、mc×mc对角阻尼比矩阵和mc×mc单位矩阵。它们随刀具进给而变化。需要注意的是，阻尼比矩阵ζc,可使用直接模态阻尼或比例阻尼建模，通常通过实验进行识别。
如上所述，刀具-主轴系统的动力学参数可以通过实验模态分析获得，而工件动力学参数是时变的，因为材料去除和刀具位置的变化会影响由ωw和Uw表示的工件动力学参数。为了准确预测铣削过程中的颤振，必须在沿刀具路径的nCS个位置对IPW的动力学参数进行采样。如果需要预测大量切削位置的IPW动力学参数，则必须对有限元模型进行耗时的网格重划分和模态再分析。
实施例2
本实施例开发了一个用于颤振预测的变IPW动力学参数的替代模型。通过将高斯过程回归与适当的正交分解相结合，确定了包括固有频率和振型在内的IPW动力学参数。第四节介绍了该方法的准确性。该方法被集成到薄壁工件铣削过程的动态模型中。实验测试和现有方法的预测结果验证了所提出的方法和动力学参数模型。最后得出结论。该方法的计算时间主要由预测一组已知的IPW动力学参数和训练GPR模型的时间组成，而预测新的IPW动力学参数花费的时间很少。通过这种方式，可以提高计算效率。
在铣削过程的动力学参数模型中，IPW动力学参数包括mw阶模态。第k阶模态由一个固有频率ω(k)和一个模态振型u(k)组成,其模态取决于刀具位置和前一刀具路径的高度，即向量e，这些映射可以表示如下：
F:e→ω(k) (5)
G:e→u(k) (6)
要准确地预测薄壁件铣削过程的稳定性叶瓣图，关键在于准确有效地识别该映射。为此，现有的方法是基于有限元计算的，而本发明将采用机器学习方法来建立一个替代模型。首先，利用现有的方法可以获得某些切削位置处e的IPW动力学参数，然后，使用高斯过程回归方法从上述已知集合中学习等式(5)中的映射。由于模态振型u(k)是一个具有大量条目的的向量，高斯过程回归不能直接用来学习映射。在高斯过程回归之前，利用POD方法对模态振型u(k)进行预处理，详细信息如下。
3.1基于高斯过程回归的变固有频率预测模型
高斯过程回归(GPR)模型是基于非参数核的概率模型。考虑已知的与切削位置对应的固有频率集合
Figure BDA0003500594530000101
其中ei＝[xi,zi,z′i]∈R3,
Figure BDA0003500594530000102
来自一个未知的分布。GPR模型解决了预测响应变量的问题即新的固有频率
Figure BDA0003500594530000103
给定新的输入向量如切削位置e*，以及训练数据，线性回归模型如下:
ω(k)＝eTβ+ε (7)
其中ε遵循一个独立的、均匀分布的高斯分布，其均值为零方差为σ2，即
Figure BDA0003500594530000104
根据数据估计了误差方差σ2和系数β。GPR模型通过从高斯过程(GP)中引入潜在变量h(ei),i＝1,2,...,nk和显式基函数φ来解释响应。潜在变量的协方差函数获取响应的平滑度，基函数将输入的e投影到p维特征空间中。
高斯过程是一组随机变量，因此它们中的任何有限数量都有一个联合高斯分布。如果{h(e),e∈R3}是一个高斯过程，则给定nk个切削位置
Figure BDA0003500594530000105
随机变量
Figure BDA0003500594530000111
的联合分布是高斯分布。高斯过程由其均值函数m(e)和协方差函数κ(e,e′)定义，即如果{h(e),e∈R3}是一个高斯过程，便可得到E(h(e))＝m(e)和Cov(h(e),h(e′))＝E[{h(e)-m(e)},{h(e′)-m(e′)}]＝κe,e′)。
现考虑已下模型：
φ(e)Tβ+h(e), (8)
其中
Figure BDA0003500594530000112
来自于一个协方差函数为κ(e,e′)的零均值高斯过程。φ(e)是一组基函数，它将R3中的切削位置向量e转换为Rp中的一个新的特征向量φ(e)。β表示基函数的系数为p维向量。此模型代表了一个高斯过程回归模型。响应的实例ω(k)可以建模为
Figure BDA0003500594530000113
因此，GPR模型是一个概率模型，对于每个切削位置ei都引入一个潜在变量h(ei)，使得GPR模型非参数化。在向量形式下，这个模型等价于
Figure BDA0003500594530000114
其中
Figure BDA0003500594530000115
潜在变量
Figure BDA0003500594530000116
在GPR模型中的联合分布如下
Figure BDA0003500594530000117
接近于线性回归模型，其中K(E,E)为
Figure BDA0003500594530000118
协方差函数κ(e,e′)也被称为核函数，它指定了两个潜在变量h(ei)和h(ei′)之间的协方差。它可以由各种函数定义，如平方指数核、指数核、有理二次核或Matren类核函数，它通常由一组核参数或超参数θ进行参数化。通常κ(e,e′)被写成κ(e,e′∣θ)来明确表示对的θ依赖性。协方差函数和基函数的形式和参数可以通过使用现有方法的超参数优化来选择。
基于特征正交分解和高斯过程回归的振型预处理模型进行变模态振型的预测：
如上所述，在高斯过程回归之前，必须使用POD方法对模态振型进行预处理。本小节将详细介绍这种预处理方法。对于特定的切削位置ei，工件的第k阶模态的模态振型表示为
Figure BDA0003500594530000121
这组模态振型的平均值为
Figure BDA0003500594530000122
通过关注与平均值的偏差，形成一组新向量：
Figure BDA0003500594530000123
假设已找到一组基函数
Figure BDA0003500594530000124
系统中任何已知的模态振型都可以精确地表示为基函数的线性组合，如下所示：
Figure BDA0003500594530000125
其中
Figure BDA0003500594530000126
最后，需要预测的任何模态振型
Figure BDA0003500594530000127
都可以近似地表示为：
Figure BDA0003500594530000128
其中，
Figure BDA0003500594530000129
是保留的基函数的数量。
关键问题是如何找到基函数
Figure BDA00035005945300001210
保留的基函数个数
Figure BDA00035005945300001211
和系数
Figure BDA00035005945300001212
设矩阵W(k)表示这些模态振型
Figure BDA00035005945300001213
的集合：
Figure BDA0003500594530000131
基函数
Figure BDA0003500594530000132
可以通过矩阵W(k)的奇异值分解得到：
W(k)＝L(k)Σ(k)R(k)T (17)
其中，
Figure BDA0003500594530000133
为左奇异矩阵。Σ(k)是奇异值矩阵，可以表示为
Figure BDA0003500594530000134
其中，
Figure BDA0003500594530000135
为第j个左奇异向量
Figure BDA0003500594530000136
的动能。
让
Figure BDA0003500594530000137
如果
Figure BDA0003500594530000138
δ是一个很小的值，如10-3。这表明，L(k)中的第一个左奇异向量
Figure BDA0003500594530000139
获取了最大可能的动能。因此，确定了基函数
Figure BDA00035005945300001310
和保留的基函数的数量
Figure BDA00035005945300001311
经过上述基函数
Figure BDA00035005945300001312
和保留基函数的个数
Figure BDA00035005945300001313
的确定过程，公式(13)中的系数
Figure BDA00035005945300001314
可以用公式(14)来计算。因此得到了一个新的映射：
Figure BDA00035005945300001315
相比等式(6)中的u(k)，α(k)显著减少。利用高斯过程回归对每个
Figure BDA0003500594530000141
的新映射(21)进行回归是可行的。
一旦得到了
Figure BDA0003500594530000142
的高斯过程回归模型，就可以计算出与新的切削位置e*相关的系数
Figure BDA0003500594530000143
最后，可以通过等式(15)计算出相应的模态振型。
为了验证所提方法的有效性和准确性，本节中将其与现有方法进行比较。主要研究了一个曲面工件，其尺寸如图2所示。曲面壁板是通过沿法线方向和相反方向偏移所需的三次贝塞尔曲线生成的。三次贝塞尔曲面描述如下：
Figure BDA0003500594530000144
工件材料为铝合金6061-T6，杨氏模量为68.9GPa，密度为2700kg/m3，泊松比为0.33。使用有限元软件对工件进行网格划分，并生成单元质量和刚度矩阵。工件的底部施加约束。网格单元类型选择为六面体单元。仿真中采用了自由度分别为31050、61440和145800的三种有限元模型。
基于现有方法的有限元模型用于计算5030个切削位置的IPW动力学参数。在所有这些模型中，刀具轨迹分为19段，以获得20个刀具位置。切削位置的z坐标随有限元z方向尺寸的增加而变化，即1.95mm。IPW动力学参数预测的计算代码用MATLAB编写，并在台式计算机上执行(Intel Xeon W-2123处理器，3.60GHz，32GB)。使用不同单元尺寸进行IPW动力学参数预测的计算时间分别为7636.12s、66510.61和458258.06，如图3所示。
所提方法采用相同的有限元模型预测IPW动力学参数。三个变量x、z和z'中的采样点数量减少一半。首先使用现有方法预测了726个切削位置的IPW动力学参数，使用不同单元尺寸的计算时间分别为1152.48s、10485.56s和59926.86s。然后，基于这些已知的IPW动力学参数，直接利用Matlab回归学习训练固有频率和模态振型的高斯过程回归模型。选择基函数和核函数作为常数函数和Matern核，参数为5/2，每个预测器分别有一个单独的长度尺度。对于振型模型，等式(20)中的δ设置为10-3。对于不同自由度的有限元模型，模态振型预处理和GPR模型训练的计算时间为1064.41s、1209.44s和1035.85s。对于三种不同自由度的有限元模型，该程序的计算时间相差不大。一旦建立了不同模态的模型，就可以使用它们来预测1ms内任意其他切削位置的IPW动力学参数。这意味着该方法的计算时间主要由预测IPW动力学参数训练集的时间和训练GPR模型的时间组成，而预测新的IPW动力学参数只需花费很少的时间。因此，对于该测试工件，与现有方法相比，该方法至少可以节省71.0％的计算时间。此外，对于具有大量自由度的有限元模型，总计算时间取决于预测已知IPW动力学参数集的计算时间。然后，为了进一步提高所提方法的计算效率，我们可以采用适用于大型薄壁件的方法。
GPR模型用于预测相同的5030个切削位置的IPW动力学参数。将前三阶固有频率与现有方法预测的频率进行了比较，如图4。该方法预测的固有频率与现有方法预测的固有频率相对误差小于3％。同时，两种方法预测的振型之间的相关性也如图5所示。主要计算两个振型向量之间夹角的余弦。三个模态对应的结果接近于1，这表明两个模态振型向量之间的夹角接近于0。这些对比意味着所提方法预测的IPW动力学参数与现有方法预测的IPW动力学参数具有可比性。
为了验证所提方法对颤振预测的有效性，使用预测的IPW动力学参数进一步预测铣削过程的稳定性。通过正斜切变换估算切削力系数，其中6061-T6的正交参数参考现有方法。结果表明，切向切削力系数为972.0N/mm2，径向为307.3N/mm2。与前三阶模态相关的工件阻尼比分别为4.99、3.97和2.46，假设它们在铣削过程中是恒定的，采用扩展的半离散时域方法对SLD进行了预测。由于在铣削过程中，工件动力学参数沿刀具路径变化，因此在经典SLD中添加了“刀具位置”作为附加尺寸，以根据刀具位置、主轴速度和测试工件的轴向切深构建三维SLD。预测的三维SLD结果如图6所示。
通过将所提方法得到的IPW动力学参数预测SLD与现有方法的实验和预测结果进行比较，验证了该方法的准确性。采用了现有方法的刀具动力学参数，预测和实验稳定性结果的比较如图7所示。注意，它们是图6的部分。在图7中，预测结果和测量结果在特定主轴转速下进行了比较。实验中，通过测量声音的频谱分布检测颤振。如果发生颤振，频谱应在刀具或工件的固有频率附近具有峰值，但与主轴或刀齿通过频率及其谐波不同。当声音或位移信号的频率分布在主轴或刀齿通过频率及其谐波上时，切削过程是稳定的。预测的稳定性极限与实验结果吻合良好。在图7所示的126组详细比较中，只有10组预测与测量结果不一致，一致率超过92％。这意味着所提出的方法可以准确地预测IPW动力学参数，并且可以集成到铣削薄壁工件的动力学参数模型中。
以上公开的本发明优选实施例只是用于帮助阐述本发明。优选实施例并没有详尽叙述所有的细节，也不限制该发明仅为所述的具体实施方式。显然，根据本说明书的内容，可作很多的修改和变化。本说明书选取并具体描述这些实施例，是为了更好地解释本发明的原理和实际应用，从而使所属技术领域技术人员能很好地理解和利用本发明。本发明仅受权利要求书及其全部范围和等效物的限制。

Claims (3)

1.一种基于高斯过程回归的薄壁结构切削颤振预测方法，其特征在于，该方法包括以下步骤：
首先，利用现有的方法获得某些切削位置向量e处的加工中工件(IPW)动力学参数，其中切削位置向量e包括当前刀心点位置参数[x，z]T和前一刀具路径高度z′，表示为e＝[x，z，z′]T，x和z分别为当前刀心点的X方向和Z方向坐标值；
然后，利用特征正交分解(POD)方法对加工中工件(IPW)动力学参数中的模态振型u(k)进行预处理，使用高斯过程回归方法从已知集合中学习映射F：e→ω(k)中的映射关系，ω(k)为固有频率；
根据材料去除和切削接触位置变化而引起的薄壁件动力学参数特性的变化，建立薄壁件铣削的动力学参数模型；所述动力学参数模型为结合高斯过程回归(GPR)和特征正交分解(POD)方法用于预测加工中工件(IPW)时变动力学参数的模型，是基于非参数核的概率模型；考虑已知的与切削位置对应的固有频率集合
Figure FDA0003500594520000011
其中
Figure FDA0003500594520000012
来自一个未知的分布；
所述动力学参数模型基于机器学习方法建立，所述机器学习方法包括：给定新的输入向量e*，以及训练数据，线性回归模型如下：
ω(k)＝eTβ+ε；
其中ε遵循一个独立的、均匀分布的高斯分布，其均值为零方差为σ2，即
Figure FDA0003500594520000013
根据数据估计了误差方差σ2和系数β；高斯过程回归(GPR)模型通过从高斯过程(GP)中引入潜在变量h(ei)，i＝1，2，...，nk和显式基函数φ来解释响应；潜在变量的协方差函数获取响应的平滑度，基函数将输入的切削位置向量e投影到p维特征空间中；
高斯过程是一组随机变量，因此它们中的任何有限数量都有一个联合高斯分布；如果{h(e)，e∈R3}是一个高斯过程，则给定nk个切削位置
Figure FDA0003500594520000014
随机变量
Figure FDA0003500594520000015
的联合分布是高斯分布；高斯过程由其均值函数m(e)和协方差函数κ(e，e′)定义，即如果{h(e)，e∈R3}是一个高斯过程，便可得到E(h(e))＝m(e)和Coy(h(e)，h(e′))＝E[{h(e)-m(e)}，{h(e′)-m(e′)}]＝κ(e，e′)；
现考虑以下模型：
φ(e)Tβ+h(e)，
其中
Figure FDA0003500594520000021
h(e)来自于一个协方差函数为κ(e，e′)的零均值高斯过程；φ(e)是一组基函数，它将变量空间R3中的切削位置向量e转换为新变量空间Rp中的一个新的特征向量φ(e)；β表示基函数的系数为p维向量；此模型代表了一个高斯过程回归模型，响应的实例ω(k)建模为：
Figure FDA0003500594520000022
高斯过程回归(GPR)模型是一个概率模型，对于每个切削位置ei都引入一个潜在变量h(ei)，使得高斯过程回归(GPR)模型非参数化；在向量形式下，这个高斯过程回归(GPR)模型等价于
Figure FDA0003500594520000023
其中
Figure FDA0003500594520000024
潜在变量
Figure FDA0003500594520000027
在高斯过程回归(GPR)模型中的联合分布如下：
Figure FDA0003500594520000025
其中K(E，E)为
Figure FDA0003500594520000026
协方差函数κ(e，e′)也被称为核函数，它指定了两个潜在变量h(ei)和h(ei′)之间的协方差；基于特征正交分解和高斯过程回归的振型预处理模型进行变模态振型的预测。
2.根据权利要求1所述的一种基于高斯过程回归的薄壁结构切削颤振预测方法，其特征在于，所述特征正交分解(POD)方法包括：
对于特定的切削位置ei，工件的第k阶模态的模态振型表示为
Figure FDA0003500594520000031
(i＝1，2，...，nk)；这组模态振型的平均值为
Figure FDA0003500594520000032
通过关注第k阶模态的模态振型与平均值的偏差，形成一组新向量：
Figure FDA0003500594520000033
假设已找到一组基函数
Figure FDA0003500594520000034
系统中任何已知的模态振型都能够精确地表示为基函数的线性组合，如下所示：
Figure FDA0003500594520000035
其中
Figure FDA0003500594520000036
最后，需要预测的任何模态振型
Figure FDA0003500594520000037
都可以近似地表示为：
Figure FDA0003500594520000038
其中，
Figure FDA0003500594520000039
是保留的基函数的数量，
Figure FDA00035005945200000310
为基函数，
Figure FDA00035005945200000311
为保留的基函数个数，
Figure FDA00035005945200000312
为系数；设矩阵W(k)表示这些模态振型
Figure FDA00035005945200000313
的集合：
Figure FDA00035005945200000314
基函数
Figure FDA0003500594520000041
可以通过矩阵W(k)的奇异值分解得到：
W(k)＝L(k)∑(k)R(k)T，
其中，
Figure FDA0003500594520000042
为左奇异矩阵，∑(k)是奇异值矩阵，可以表示为
Figure FDA0003500594520000043
其中，
Figure FDA00035005945200000414
为第j个左奇异向量
Figure FDA0003500594520000044
的动能；
让
Figure FDA0003500594520000045
如果
Figure FDA0003500594520000046
δ≤10-3，L(k)中的第一个左奇异向量
Figure FDA0003500594520000047
获取了最大可能的动能，确定了基函数
Figure FDA0003500594520000048
和保留的基函数的数量
Figure FDA0003500594520000049
经过所述基函数
Figure FDA00035005945200000410
和保留基函数的个数
Figure FDA00035005945200000411
的确定过程，得到了一个新的映射：
Figure FDA00035005945200000412
计算出与新的切削位置向量e*相关的系数
Figure FDA00035005945200000413
最后，通过等式
Figure FDA0003500594520000051
计算出相应的模态振型。
3.根据权利要求1-2所述的一种基于高斯过程回归的薄壁结构切削颤振预测方法，其特征在于，所述协方差函数κ(e，e′)由平方指数核、指数核、有理二次核或Matren类核函数定义，由一组核参数或超参数θ进行参数化。
CN202210126159.2A 2022-02-10 2022-02-10 一种基于高斯过程回归的薄壁结构切削颤振预测方法 Active CN114626015B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210126159.2A CN114626015B (zh) 2022-02-10 2022-02-10 一种基于高斯过程回归的薄壁结构切削颤振预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210126159.2A CN114626015B (zh) 2022-02-10 2022-02-10 一种基于高斯过程回归的薄壁结构切削颤振预测方法

Publications (2)

Publication Number Publication Date
CN114626015A true CN114626015A (zh) 2022-06-14
CN114626015B CN114626015B (zh) 2024-03-01

Family

ID=81897969

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210126159.2A Active CN114626015B (zh) 2022-02-10 2022-02-10 一种基于高斯过程回归的薄壁结构切削颤振预测方法

Country Status (1)

Country Link
CN (1) CN114626015B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115544656A (zh) * 2022-09-30 2022-12-30 华中科技大学 一种薄壁叶片加工时变模态参数高效预测方法及系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106096146A (zh) * 2016-06-15 2016-11-09 西北工业大学 切削过程中薄壁件动力学参数的预测方法
CN107505842A (zh) * 2017-09-04 2017-12-22 重庆邮电大学 一种数控机床广义空间切削稳定性预测与优化方法
WO2020218982A1 (en) * 2019-04-24 2020-10-29 Sabanci Üniversitesi A method for generating a tool path to manufacture a part using a computer numerical control machine system

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106096146A (zh) * 2016-06-15 2016-11-09 西北工业大学 切削过程中薄壁件动力学参数的预测方法
CN107505842A (zh) * 2017-09-04 2017-12-22 重庆邮电大学 一种数控机床广义空间切削稳定性预测与优化方法
WO2020218982A1 (en) * 2019-04-24 2020-10-29 Sabanci Üniversitesi A method for generating a tool path to manufacture a part using a computer numerical control machine system

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
张雪薇;于天彪;王宛山;: "薄壁零件铣削三维颤振稳定性建模与分析", 东北大学学报(自然科学版), no. 01, 15 January 2015 (2015-01-15) *
胡耀斌;邹湘军;张春良;江涌涛;: "SVR在切削颤振状态趋势预测中的应用", 系统仿真学报, no. 17, 5 September 2007 (2007-09-05) *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115544656A (zh) * 2022-09-30 2022-12-30 华中科技大学 一种薄壁叶片加工时变模态参数高效预测方法及系统

Also Published As

Publication number Publication date
CN114626015B (zh) 2024-03-01

Similar Documents

Publication Publication Date Title
Adetoro et al. An improved prediction of stability lobes using nonlinear thin wall dynamics
Liao et al. Tool wear state recognition based on GWO–SVM with feature selection of genetic algorithm
Chen et al. Automated surface subdivision and tool path generation for 31212-axis CNC machining of sculptured parts
Ramesh et al. Support vector machines model for classification of thermal error in machine tools
Lin et al. Optimization with minimum process error for layered manufacturing fabrication
Li et al. Time-varying dynamics updating method for chatter prediction in thin-walled part milling process
Fountas et al. An integrated framework for optimizing sculptured surface CNC tool paths based on direct software object evaluation and viral intelligence
Liang et al. An accuracy algorithm for chip thickness modeling in 5-axis ball-end finish milling
Wang et al. A prediction method based on the voxel model and the finite cell method for cutting force-induced deformation in the five-axis milling process
Huang et al. A systematic approach for online minimizing volume difference of multiple chambers in machining processes based on high-definition metrology
Sun et al. Iso-planar feed vector-fields-based streamline tool path generation for five-axis compound surface machining with torus-end cutters
Su et al. An image-based approach to predict instantaneous cutting forces using convolutional neural networks in end milling operation
CN114626015B (zh) 一种基于高斯过程回归的薄壁结构切削颤振预测方法
Xi et al. A prediction model of the cutting force–induced deformation while considering the removed material impact
Lei et al. Prediction of the posture-dependent tool tip dynamics in robotic milling based on multi-task Gaussian process regressions
CN116520770A (zh) 五轴联动数控机床伺服动态特性匹配特征值评价方法
Huang et al. Collision detection algorithm on abrasive belt grinding blisk based on improved octree segmentation
Wu et al. Computer aided machining fixture design algorithm and software based on case learning for near-net-shaped jet engine blade
Van Dang et al. Enhanced vector flow of significant directions for five-axis machining of STL surfaces
Yang et al. CNC corner milling parameters optimization based on variable-fidelity metamodel and improved MOPSO regarding energy consumption
Franco et al. Optimal cutting condition selection for high quality receptance measurements by sweep milling force excitation
Wu et al. Layout optimization of auxiliary support for deflection errors suppression in end milling of flexible blade
Nikolaev et al. Optimal milling modes identification of a jet-engine blade using time-domain technique
Hu et al. Error prediction of balancing machine calibration based on machine learning method
CN113570147B (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
GR01 Patent grant
GR01 Patent grant