CN102564787A - 基于空运行激励的数控机床模态比例因子获取方法 - Google Patents

Abstract

本发明公开了一种基于空运行激励的数控机床模态比例因子获取方法，包括以下步骤：生成机床加速度的二值随机序列，并根据二值随机序列生成机床的空运行数控代码，计算二值随机惯性激励力序列的自功率谱G_XX(jω)，执行空运行数控代码，以测量机床的响应信号并计算响应信号的互功率谱矩阵，根据响应信号的互功率谱矩阵利用最小二乘复频域法计算系统极点λ_1...N和

以及模态振型向量ψ_1...N和

根据自功率谱[G_XX(jω)]、系统极点λ_1...N和

以及模态振型向量ψ_1...N和

计算机床结构的模态比例因子。本发明能够估计激励序列的能量大小，继而从机床测点间的互功率谱矩阵中获取模态比例因子。

Description

基于空运行激励的数控机床模态比例因子获取方法

技术领域

本发明涉及数控装备结构模态参数分析技术领域，尤其涉及一种基于空运行激励的数控机床模态比例因子获取方法。

背景技术

传统实验模态分析技术能够得到完整的模态参数，包括固有频率、阻尼比、模态振型向量和模态比例因子。但对于大型结构(如桥梁、高楼和重型机床)，传统实验模态分析技术的激励方式难以实施，利用其它易于实施的随机激励(如环境激励)作为激励方式的工作模态分析方法得到了重视。工作模态分析方法中，假设输入激励为幅值恒定的白噪声随机激励，推导出测点间的互功率谱矩阵与频响函数矩阵有相似的表达式，可通过互功率谱矩阵辨识出部分的模态参数。

现有的工作模态分析方法中，能够辨识结构的固有频率、阻尼比和模态振型向量这些模态参数，但无法得到模态比例因子。由于现有工作模态激励属于环境激励，无法控制激励输入的位置，并且难以测量和估计实际激励力的能量大小，从而无法直接从互功率谱矩阵中得到模态比例因子，不能得到结构完整的频响函数。

发明内容

本发明的目的在于提供一种基于空运行激励的数控机床模态比例因子获取方法，其能够估计激励序列的能量大小，继而从机床测点间的互功率谱矩阵中获取模态比例因子。

本发明是通过以下技术方案实现的：

一种基于空运行激励的数控机床模态比例因子获取方法，包括以下步骤：

(1)生成机床加速度的二值随机序列，并根据二值随机序列生成机床的空运行数控代码：

(2)计算二值随机惯性激励力序列的自功率谱G_XX(jω)：

(3)执行空运行数控代码，以测量机床的响应信号并计算响应信号的互功率谱矩阵；

(4)根据响应信号的互功率谱矩阵利用最小二乘复频域法计算系统极点λ_1...N和

以及模态振型向量ψ_1...N和

(5)根据自功率谱[G_XX(jω)]、系统极点λ_1...N和

以及模态振型向量ψ_1...N和计算机床结构的模态比例因子。

步骤(1)包括以下子步骤：

(2-1)根据机床模态分析所需的频带确定钟频率f₀的大小，f₀为频带的两倍；

(2-2)确定机床在空运行过程中产生的惯性力的绝对值ma，其中a为机床的加速度，m为机床的质量；

(2-3)通过MATLAB软件生成钟频率为f₀、二值变量在0和ma中取值的二值随机惯性激励力序列；

(2-4)根据二值随机惯性激励力序列计算生成加速度二值随机序列。

(2-5)根据加速度二值随机序列编写机床的空运行数控代码；

(2-6)根据二值随机序列生成机床的空运行数控代码。

步骤(2)是采用以下公式：

$G_{\mathrm{xx}}\left(\mathrm{j\ω}\right)=\frac{{\left(\mathrm{ma}\right)}^2}{2f_0}{\left(\frac{\sin \mathrm{\ω}/f_0}{\mathrm{\ω}/f_0}\right)}^2,$

其中G_xx(jω)为二值随机惯性激励力序列的自功率谱函数。步骤(3)是采用以下公式：

$\left[G_{\mathrm{yy}}\left(\mathrm{j\ω}\right)\right]=\left[\begin{array}{cccc}{\mathrm{PSD}}_{11}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{12}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{13}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{14}\left(\mathrm{j\ω}\right)\\ {\mathrm{CSD}}_{21}\left(\mathrm{j\ω}\right)& {\mathrm{PSD}}_{22}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{23}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{24}\left(\mathrm{j\ω}\right)\\ {\mathrm{CSD}}_{31}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{32}\left(\mathrm{j\ω}\right)& {\mathrm{PSD}}_{33}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{34}\left(\mathrm{j\ω}\right)\\ {\mathrm{CSD}}_{41}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{42}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{43}\left(\mathrm{j\ω}\right)& {\mathrm{PSD}}_{44}\left(\mathrm{j\ω}\right)\end{array}\right]---\left(3\right)$

其中[G_yy(jω)]表示响应信号的互功率谱矩阵，CSD_kl(jω)表示第k点响应信号对l点响应信号的互功率谱，PSD_kk(jω)表示第k点响应信号的自功率谱。

本发明具有以下的优点和技术效果：

本发明所述的激励只需设定机床的加速度并生成数控机床的空运行数控代码以使机床工作台随机往复运动，该往复运动可以产生二值随机惯性激励力序列，利用二值随机序列的频谱公式计算出该随机输入的能量大小，然后通过随机输入的能量大小计算获得频响函数的模态比例因子，从而能够得到整个机床完整的频响函数。

附图说明

图1为本发明基于空运行激励的数控机床模态比例因子获取方法的流程图。

图2为本发明二值随机惯性激励力序列的示意图。

图3为本发明二值随机惯性激励力序列的时域图和自功率频谱图。

具体实施方式

首先对本发明的技术术语进行解释和说明：

二值随机惯性激励力序列：机床工作台加减速过程中产生的惯性力按照二值随机序列分布形成的激励序列。

空运行数控代码：用于控制数控机床实现空运行的指令代码。

以下以XHK5140型机床为例对本发明进行说明。

如图1所示，本发明基于空运行激励的数控机床模态比例因子获取方法包括以下步骤：

(1)生成机床加速度的二值随机序列，并根据二值随机序列生成机床的空运行数控代码：

具体而言，根据二值随机序列的定义，二值随机序列包含两个参数：钟频率f₀和二值变量的取值范围。二值随机序列的自功率谱如图2所示，在0-0.5f₀范围内基本保持平直，可将其视为“白噪声”，Δt表示序列的单位时间间隔，Δt＝1/f₀，φ表示初始相位。

根据机床模态分析所需的频带确定钟频率f₀的大小，f₀为机床模态分析所需带宽的两倍。

设置机床的加速度控制参数，使机床工作台空运行的加速度值最大，该加速度的绝对值为a，则机床工作台(质量为m)在空运行过程中产生的惯性力的绝对值为ma。

通过MATLAB软件生成钟频率为f₀、二值变量在0和ma中取值的二值随机惯性激励力序列，该序列的示意图见图3，图中n₁、n₂为大于或等于零的正整数。

根据二值随机惯性激励力序列计算生成加速度二值随机序列，然后利用速度、加速度、位移和时间之间的计算公式，编写机床的空运行数控代码。

(2)计算二值随机惯性激励力序列的自功率谱G_XX(jω)：

具体而言，根据二值随机序列的自功率谱计算公式有：

$G=\frac{{2A}^2}{f_0}{\left(\frac{\sin \mathrm{\ω}/f_0}{\mathrm{\ω}/f_0}\right)}^2---\left(1\right)$

其中，G为二值随机序列的自功率谱函数，A为二值随机序列变量最大值与最小值之差的绝对值的一半。该频谱在0-0.5f₀范围内基本保持平直，可将其视为“白噪声”，其能量大小为：

将对应的二值随机惯性激励力序列的参数代入公式(1)中，可以计算出输入激励的自功率谱为：

$G_{\mathrm{xx}}\left(\mathrm{j\ω}\right)=\frac{{\left(\mathrm{ma}\right)}^2}{2f_0}{\left(\frac{\sin \mathrm{\ω}/f_0}{\mathrm{\ω}/f_0}\right)}^2---\left(2\right)$

其中G_xx(jω)为输入的二值随机惯性激励力序列的自功率谱函数。该频谱在0-0.5f₀范围内基本保持平直，可将其视为“白噪声”，其能量大小为：

(3)执行空运行数控代码，以测量机床的响应信号并计算响应信号的互功率谱矩阵；

具体而言，根据以下公式计算响应信号的互功率谱矩阵(以机床的4个响应点为例)：

$\left[G_{\mathrm{yy}}\left(\mathrm{j\ω}\right)\right]=\left[\begin{array}{cccc}{\mathrm{PSD}}_{11}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{12}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{13}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{14}\left(\mathrm{j\ω}\right)\\ {\mathrm{CSD}}_{21}\left(\mathrm{j\ω}\right)& {\mathrm{PSD}}_{22}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{23}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{24}\left(\mathrm{j\ω}\right)\\ {\mathrm{CSD}}_{31}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{32}\left(\mathrm{j\ω}\right)& {\mathrm{PSD}}_{33}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{34}\left(\mathrm{j\ω}\right)\\ {\mathrm{CSD}}_{41}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{42}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{43}\left(\mathrm{j\ω}\right)& {\mathrm{PSD}}_{44}\left(\mathrm{j\ω}\right)\end{array}\right]---\left(3\right)$

其中[G_yy(jω)]表示响应信号的互功率谱矩阵，GSD_kl(jω)表示第k点响应信号对l点响应信号的互功率谱，PSD_kk(jω)表示第k点响应信号的自功率谱。

(4)根据响应信号的互功率谱矩阵利用最小二乘复频域法计算系统极点λ_1...N和以及模态振型向量ψ_1...N和

根据结构频响函数的特性可知：

[G_yy(jω)]＝[H(jω)]^*[G_XX(jω)][H(jω)]^T          (4)

其中[G_yy(jω)]表示响应信号的互功率谱矩阵，[H(jω)]表示频响函数矩阵，[G_XX(jω)]表示激励力的自功率谱。

工作模态分析方法中假设激励力的自功率谱为白噪声，则G_XX(jω)＝常数C，根据频响函数的留数形式：

$\left[H\left(\mathrm{j\ω}\right)\right]=\underset{r=1}{\overset{N}{\mathrm{\Σ}}}(\frac{Q_r{\mathrm{\ψ}}_r{\mathrm{\ψ}}_r^T}{\mathrm{j\ω}-{\mathrm{\λ}}_r}+\frac{Q_r^{*}{\mathrm{\ψ}}_r^{*}{\mathrm{\ψ}}_r^{*T}}{\mathrm{j\ω}-{\mathrm{\λ}}_r^{*}})---\left(5\right)$

$=\left[\mathrm{\Ψ}\right]\left[\mathrm{j\ω}\right[I]-[\mathrm{\Λ}\left]\right]\left[Q\right]{\left[\mathrm{\Ψ}\right]}^T$

其中，[H(jω)]表示频响函数矩阵，

ψ_1...N和

是模态振型向量，[Λ]是以系统极点λ_1...N和

组成的对角阵，[Q]是以频响函数的模态比例因子Q_1...N和

组成的对角阵。

将公式(5)和G_XX(jω)＝C代入公式(4)得

[G_yy(jω)]＝[Ψ]^*[jω[I]-[Λ]]^*[Q]^*[Ψ]^HC[Ψ][Q]^T[jω[I]-[Λ]]^T[Ψ]^T

                                                            (6)

若设

[D]＝[Q]^*[Ψ]^HC[Ψ][Q]^T[jω[I]-[Λ]]^T                      (7)

由于C为常数，根据模态振型向量的定义，[Ψ]^HC[Ψ]仍为对角阵，而且[jω[I]-[Λ]]^T也是一个对角阵，因此[D]也保持了对角阵的形式。设[D]为以

和d_1...N组成的对角阵。

则

$\left[G_{\mathrm{yy}}\left(\mathrm{j\ω}\right)\right]={\left[\mathrm{\Ψ}\right]}^{*}\left[\mathrm{j\ω}\right[I]-[\mathrm{\Λ}\left]{]}^{*}\right[D]{\left[\mathrm{\Ψ}\right]}^T$

$=\underset{r=1}{\overset{N}{\mathrm{\Σ}}}(\frac{d_r{\mathrm{\ψ}}_r{\mathrm{\ψ}}_r^T}{\mathrm{j\ω}-{\mathrm{\λ}}_r}+\frac{d_r^{*}{\mathrm{\ψ}}_r^{*}{\mathrm{\ψ}}_r^{*T}}{\mathrm{j\ω}-{\mathrm{\λ}}_r^{*}})---\left(8\right)$

此时，可见[G_yy(jω)和[H(jω)]有相同的表达形式，因此用[G_yy(jω)代替[H(jω)]用最小二乘复频域法(LSCE)可以得到公式(8)中的λ_1...N，

ψ_1...N，

d_1...N，

其中d_1...N和

为互功率谱的模态比例因子；

(5)根据前述计算得到的二值随机惯性激励力序列的自功率谱[G_XX(jω)]、系统极点λ_1...N和

以及模态振型向量ψ_1...N和

计算机床结构的模态比例因子：

利用步骤(2)里计算得到的二值随机惯性激励力序列的自功率谱[G_XX(jω)]，然后将常数C和步骤(4)中计算得到的[D]、[Ψ]和[Λ]代入公式(7)，可解得频响函数的模态比例因子Q_1...N和

Claims (4)

1.一种基于空运行激励的数控机床模态比例因子获取方法，其特征在于，包括以下步骤：

(1)生成机床加速度的二值随机序列，并根据所述二值随机序列生成机床的空运行数控代码：

(2)计算所述二值随机惯性激励力序列的自功率谱G_XX(jω)：

(3)执行所述空运行数控代码，以测量机床的响应信号并计算所述响应信号的互功率谱矩阵；

(4)根据所述响应信号的互功率谱矩阵利用最小二乘复频域法计算系统极点λ_1...N和

以及模态振型向量ψ_1...N和

(5)根据所述自功率谱[G_XX(jω)]、系统极点λ_1...N和以及模态振型向量ψ_1...N和

计算机床结构的模态比例因子。

2.根据权利要求1所述的方法，其特征在于，所述步骤(1)包括以下子步骤：

(2-1)根据机床模态分析所需的频带确定钟频率f₀的大小，f₀为所述频带的两倍；

(2-2)确定机床在空运行过程中产生的惯性力的绝对值ma，其中a为机床的加速度，m为机床的质量；

(2-3)通过MATLAB软件生成钟频率为f₀、二值变量在0和ma中取值的二值随机惯性激励力序列；

(2-4)根据所述二值随机惯性激励力序列计算生成加速度二值随机序列。

(2-5)根据所述加速度二值随机序列编写机床的空运行数控代码；

(2-6)根据所述二值随机序列生成机床的空运行数控代码。

3.根据权利要求1所述的方法，其特征在于，所述步骤(2)是采用以下公式：

$G_{\mathrm{xx}}\left(\mathrm{j\ω}\right)=\frac{{\left(\mathrm{ma}\right)}^2}{2f_0}{\left(\frac{\sin \mathrm{\ω}/f_0}{\mathrm{\ω}/f_0}\right)}^2,$

其中G_xx(jω)为所述二值随机惯性激励力序列的自功率谱函数。

4.根据权利要求3所述的方法，其特征在于，所述步骤(3)是采用以下公式：

$\left[G_{\mathrm{yy}}\left(\mathrm{j\ω}\right)\right]=\left[\begin{array}{cccc}{\mathrm{PSD}}_{11}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{12}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{13}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{14}\left(\mathrm{j\ω}\right)\\ {\mathrm{CSD}}_{21}\left(\mathrm{j\ω}\right)& {\mathrm{PSD}}_{22}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{23}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{24}\left(\mathrm{j\ω}\right)\\ {\mathrm{CSD}}_{31}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{32}\left(\mathrm{j\ω}\right)& {\mathrm{PSD}}_{33}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{34}\left(\mathrm{j\ω}\right)\\ {\mathrm{CSD}}_{41}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{42}\left(\mathrm{j\ω}\right)& {\mathrm{CSD}}_{43}\left(\mathrm{j\ω}\right)& {\mathrm{PSD}}_{44}\left(\mathrm{j\ω}\right)\end{array}\right]---\left(3\right)$

其中[G_yy(jω)]表示所述响应信号的互功率谱矩阵，CSD_kl(jω)表示第k点响应信号对l点响应信号的互功率谱，PSD_kk(jω)表示第k点响应信号的自功率谱。

Images