机床是由许多零部件按照一定要求结合起来的，零部件之间相互结合的部位称为结合部^[1]。结合部动力学特性对整机的固有频率、振型等动力学特性影响显著。研究表明，整机设备中结合部的接触刚度约占总刚度的60%～80%，结合部的接触阻尼占全部阻尼的90%以上，结合部引起的变形量约占整机总变形量的40%～60%^[2-3]。而螺栓结合部作为机床装配中的三大典型结合部之一，其动力学参数的准确辨识是整机有限元动力学建模与性能优化的前提，因此具有重要的研究价值。

为寻找比较理想的螺栓结合部动力学参数的识别方法，国内外学者做了大量研究工作，通常采用的方法有频响函数法^[4-6]和模态实验法^[7-11]。前者能有效的辨识整个激振频率下的结合部刚度和阻尼，但该方法会因辨识过程中频响函数的逆运算造成辨识结果较大误差。后者通过理论计算不断缩小理论模型与实际模型间的误差直至达到收敛，但辨识过程中调用有限元分析次数较多，辨识精度虽高但工作效率低。

目前，响应面法已在一些工程优化领域广泛的应用。响应面法是基于统计学和综合实验技术来解决复杂系统输入与输出之间映射关系的一种近似方法，即少数观测样本的统计回归，采用显式函数形式表达复杂隐式关系^[12-13]。作者基于响应面法将研究对象结合面间的刚度与动态性能指标—模态固有频率关系显现化，从而实现结合面间的参数辨识。通过响应函数的选型对结合部刚度间的耦合关系显示化，进而分析结合部刚度间耦合对组合件动态性能的影响。

1 螺栓结合部动力学模型

在对螺栓结合部模型进行模态分析时，组合件在各阶固有频率下产生共振，此时结合面刚度为动刚度。结合部动力学特性用弹簧–阻尼单元模拟，等效螺栓结合部的动力学模型如图1所示。A和B分别为螺栓结构中的子结构，K为结合部等效刚度，C为结合部等效阻尼。

螺栓结合部动力学方程为：

$ \left( \!\!\!{\begin{array}{*{20}{c}} {{{{M}}_{\rm A}}}\!\!\!\!&\!\!\!\!0\\ 0\!\!\!\!&\!\!\!\!{{{{M}}_{\rm B}}} \end{array}}\!\!\! \right)\!{{\ddot X}}\left( t \right) \!+\! \left(\!\!\! {\begin{array}{*{20}{c}} {{{{C}}_{\rm A}}}\!\!\!\!&\!\!\!\!{ - {{{C}}_{\rm J}}}\\ { - {{{C}}_{\rm J}}}\!\!\!\!&\!\!\!\!{{{{C}}_{\rm B}}} \end{array}} \!\!\!\right)\!{{\dot X}}\left( t \right) \!+\! \left(\!\!\! {\begin{array}{*{20}{c}} {{{{K}}_{\rm A}}}\!\!\!\!&\!\!\!\!{ - {{{K}}_{\rm J}}}\\ { - {{{K}}_{\rm J}}}\!\!\!\!&\!\!\!\!{{{{K}}_{\rm B}}} \end{array}}\!\!\! \right)\!{{X}}\left( t \right)\! = \!{{F}}\left( t \right) $

(1)

为方便得到结合部间的等效刚度与阻尼，将自由状态下的螺栓结合部动态特性作为研究对象，即：K_A、K_B、C_A、C_N为0。

文献[14-15]研究表明，如装配条件已定时，其结构材料、几何尺寸和装配位置是确定的，那么装配后的结构件结合部动态特性的模态固有频率ω与刚度矩阵K是一一对应的映射关系，关系如下：

${{\omega }} = function({{ k}_1},{{ k}_2}, \cdots, {{ k}_n})$

(2)

式中， $\omega $ 为模态固有频率矢量， ${{k}}$ 为结合面间各节点间的刚度矩阵。式（2）中螺栓组合件模态固有频率与其结合部刚度矩阵关系是隐性的，为解决该问题，作者基于响应面法将二者关系显示化，通过拟合其响应面函数即可预测研究对象的结合部等效动刚度参数。

2 基于响应面法的结合部参数辨识 2.1 实验设计方案

实验设计方案的选择是响应面法构建近似响应函数的前提与保证，因中心复合设计（CCD）实验设计方法是中心组合设计的一种类型，可以保持对高阶响应面定义的基础上去掉部分实验，属于简单和采样均匀分布的组合设计，本文基于参数变量的数目采用CCD方法制定响应面法的实验样本点采样。

2.2 响应面法的原理

为研究结合部刚度间的耦合关系，选择响应面函数为含交叉项二次多项式，函数模型如式（3）所示：

$\widetilde Y(x) = \sum\limits_{i = 0}^{p - 1} {{\beta _i}} {{ k}_i}$

(3)

式中， ${\beta _i}$ 为函数待定系数，因所选形式为含交叉项的二次多项式，其个数 $p$ 为 $p = (n + 1)(n + 2)/2$ ，n为结合面间不同刚度参数数目。

${{\xi }}$ 为响应面函数 ${{\widetilde Y}}$ 与真实函数 $ Y$ 的随机误差向量，可用下列矩阵表示：

${{\xi}} = \widetilde { Y}({{k) - Y}}({{k}})$

(4)

为找到最接近所有实验数据点的响应面，利用最小二乘原理使 ${{\xi }}$ 的平方和最小即满足式（5）的条件^[13]：

$\begin{aligned}[b] S({{\beta }}) &= {\sum\limits_{j = 0}^{m - 1} {({{{\xi }}^j})} ^2} = {\sum\limits_{j = 0}^{m - 1} {\left(\sum\limits_{i = 0}^{p - 1} {{\beta _i}} { k}_i^j - {y^j}\right)} ^2} \to \min, \\ \frac{{\partial S}}{{\partial {{\beta }}}} &= 2\sum\limits_{j = 0}^{m - 1} {\left[ {{ k}_i^j\left(\sum\limits_{i = 0}^{p - 1} {{\beta _i}} { k}_i^j - {y^j}\right)} \right]} = 0 \end{aligned} $

(5)

式中，m为取样次数。将式（5）化简可得回归系数向量 $\;{{\beta }}$ ：

${{\beta = (}}{{{K}}^{{T}}}{{K}}{{{)}}^{{{ - 1}}}}{{{K}}^{\rm{T}}}{{Y}}$

(6)

将式（6）代入到式（3）即可得到设计变量的近似响应函数二次多项式。

为了评价响应面函数对响应值（实验数据）拟合的程度，需要采用评价指标来评价，常用的评价指标有复相关系数 ${R^2}$ ：

${R^2} = 1{{ - }}\frac{{{{{Y}}^{\rm{T}}}{{Y - }}{{{\beta }}^{\rm{T}}}{{{X}}^{\rm{T}}}{{Y}}}}{{{{{Y}}^{\rm{T}}}{{Y - }}\dfrac{{{{({{{I}}^{\rm{T}}}{{Y)}}}^{2}}}}{m}}}$

(7)

式中， ${R^2}$ 是一个在[0,1]之间变化的值，它的值越接近1说明误差的影响越小，即响应面函数对响应值的拟合越准确。

2.3 遗传算法和非线性规划的函数寻优算法 2.3.1 目标函数与设计变量

通过M+P模态实验和响应面函数分别获得前n阶固有频率的实验值 $\omega _j^c$ 和计算值 ${\omega _j}({k_i})(j=1, 2, \cdots , {n})$ ，如果各阶 $\omega _j^c$ 与对应的各阶 ${\omega _j}({k_i})$ 差值 $\xi $ 在允许范围内，那么可以认为所选定的自变量值即为要确定的螺栓结合部的等效动刚度。于是建立目标函数如下：

$\begin{aligned} \min f(k) = \frac{1}{m}\sum\limits_{j = 1}^m \left(\frac{{{\omega _j}({k_1},{k_{\rm{2}}} ,\cdots ,{k_i}) - \omega _j^c}}{{\omega _j^c}} \right)^2 \le \xi,\\ j = 1,2, \cdots, m;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\\ \;\;\;\;\;{\rm{s.t.}}\;\;\;\;\;LB \le {k_i} \le UB,\;\;\;\;\;i = 1,2 ,\cdots, {n} \end{aligned}$

(8)

2.3.2 寻优算法选择与应用

式（8）的目标函数是关于结合部各结合点处接触刚度的非线性函数，解决此类非线性优化问题可以采用非线性算法。但此算法大多数采用梯度下降的方法求解，局部搜索能力较强，但全局搜索能力较弱。而遗传算法采用选择、交叉和编译算子进行搜索，全局搜索能力较强，但是局部搜索能力较弱，因此采用遗传算法和非线性规划的函数寻优算法可以结合两种算法的优点^[16]，遗传算法迭代计算1次，以遗传算法当前计算的结果为初始值，采用MATLAB中线性规划函数fmincon进行局部寻优，并把寻找到的局部最优值作为新个体染色体继续进化，以得到问题的全局最优解。

2.4 结合部参数辨识技术路线

基于第2.1～2.3节阐述的原理及方法可知，通过构建广义模态固有频率响应面函数，完成研究对象结合面间的刚度与各阶模态固有频率的关系分析及实现结合部刚度参数的辨识工作，其技术路线主要步骤归纳如下：

1）采用德国M+P实验设备完成螺栓组合件的自由模态动力学实验，获取前11阶模态固有频率。

2）依据子结构实验数据对子结构进行模型修正，建立准确的螺栓结构件动力学有限元分析模型，进行网格划分，并在其结合面间指定位置生成单对或多对弹簧–阻尼单元。

3）基于实验因子、因子水平个数，采用CCD实验设计方案制定实验采样点的矩阵表。

4）根据实验采样矩阵表创建多个有限元分析命令流文件，每个文件中写入每个位置的不同刚度信息，以此形成研究对象的有限元模型实验设计动力学分析样本空间。

5）有限元计算并提取计算结果，基于每组实验样本点与实验结果在MATLAB中构建响应面近似函数，并检验响应面函数的拟合度。

6）基于遗传算法和非线性规划的函数寻优算法完成结合部参数的辨识。

3 实例研究 3.1 螺栓结合部实验设计及数据采集

按照机床主轴与主轴箱的装配方式对螺栓组合件进行装配与预紧（螺栓型号：M8，预紧力：80 N·m），选用材质为45钢，子结构尺寸为mm。分别对子结构和组合结构进行模型实验，通过弹力绳悬挂模拟结构自由边界。实验设备为德国M+P振动噪声测试平台，激振方式为力锤激励法，数据采集使用的加速度传感器为PCB ICP（型号：3273A2，重量约为4 g），因与所测结构中两子结构的质量比均小于5%，所以传感器的附加质量可忽略不记。测点尽量分布于结构的中心线上，为减少随机误差的影响，每个测点的频响取值次数设置为4次，取其平均值作为最终采样频响。

利用smart office软件对所建立的几何模型及采集的所有测点的实验数据进行模态和振型分析。选择在0～3 600 Hz频率区间进行模态参数识别，数据分析后所得前11阶固有频率如表1所示。

实验配置及测点分布如图2所示。

3.2 螺栓结构件有限元模型建立

为准确建立螺栓结构件有限元模型，基于子结构的模态实验对有限元模型进行修正，修正后的前5阶总模态固有频率误差由18.1%降至4.3%，平均误差仅为0.87%。两子结构通过弹簧命令模拟螺栓结合部的动力学特性，为研究有限元模型结合面间多对节点和单对节点对组合件结合部动力学的预测，弹簧按下面两种方式布置^[17-18]：

1）法向和切向弹簧分别分布于结合面的两端中心位置，共4组弹簧，建立的螺栓结合部有限元模型如图3（a）所示。

2）法向和切向弹簧分别分布于结合面的中心位置，共2组弹簧，建立的螺栓结合部有限元模型见图3（b）。

3.3 CCD实验设计方案

按两对节点弹簧布置方式可知n=4，图3所示的参数 ${k_1}$ 、 ${k_2}$ 为法向刚度， ${k_{3}}$ 、 ${k_{4}}$ 为切向刚度。因参数的取值范围不可知，预估范围为（k₁、k₂、k₃、k₄） $ \in (1 \times {10^{3}},{\rm{1}} \times$ $ {10^{\rm{5}}}) $ N/mm^[1]。根据CCD实验设计方案，建立实验因子水平（表2）和25组仿真采样矩阵（表3）。

3.4 响应面函数的构建

对表3中每组采样点进行有限元模态分析，提取前11阶模态固有频率^[19]。在MATLAB中建立螺栓结合部刚度和广义模态固有频率函数关系式：

$\begin{aligned}[b] &{{\omega }} = {{{\beta }}_{0}} + {{{\beta }}_{1}}{k_1} + {{{\beta }}_2}{k_2} + {{{\beta }}_{3}}{k_{3}} + {{{\beta }}_{4}}{k_{4}} + {{{\beta }}_{5}}k_1^2 + {{{\beta }}_{6}}k_2^2 + {{{\beta }}_{7}}k_{3}^2 + \\ &{{{\beta }}_{8}}k_{4}^2\! +\! {{{\beta }}_{9}}{k_1}{k_2} \!+\! {{{\beta }}_{10}}{k_1}{k_{3}} \!+\! {{{\beta }}_{{11}}}{k_1}{k_{4}} \!+ \!{{{\beta }}_{{12}}}{k_{2}}k_{3}\! +\! {{{\beta }}_{{13}}}{k_{2}}{k_{4}} \!+ \!{{{\beta }}_{14}}{k_{3}}{k_{4}} \\ \end{aligned} $

(9)

式中， ${{{\beta }}_{0}} {\text{～}} {{{\beta }}_{4}}$ 为线性系数向量， ${{{\beta }}_{5}} {\text{～}}{{{\beta }}_{8}}$ 为非线性系数向量， ${{{\beta }}_{9}} {\text{～}} {{{\beta }}_{{14}}}$ 为交叉耦合系数。

计算响应面方程的15个系数矩阵如式（10）所示：

(10)

由式（10）数据可知，结合部建模的多个刚度相互耦合，最终以结合部的动刚度属性影响组合结构的动力学特性。 ${{{\beta }}_{{10}}} \sim {{{\beta }}_{{13}}}$ 系数数值较小，说明法向与切向间刚度耦合可以忽略不计； ${{{\beta }}_{9}}$ 、 ${{{\beta }}_{14}}$ 系数数值同 ${{{\beta }}_{1}} \sim {{{\beta }}_{8}}$ 数值为同量级，说明法向间和切向间耦合关系不容忽视。 ${{{\beta }}_{5}} \sim {{{\beta }}_{8}}$ 的数据说明了螺栓结合部刚度具有非线性的特性。

选取表3中的任意采样组合点共15组数据，计算复相关系数 $R^2$ ，进行多次取样测试，结果表明 $R^2$ 与1非常接近（差值均小于0.001），说明响应近似函数拟合精度已足够。

将式（10）代入到式（9）中即可得到优化后的响应面函数，通过此函数方程在参数取值范围内按照遗传算法与二次规划算法联合寻优法即可得到4刚度组合的最优值①，按上述流程对单对节点刚度进行辨识，得到最优值②。结果如表4所示。

为研究单对节点、2对节点4组刚度耦合/无耦合对有限元模型预测能力的影响，分别建立了考虑刚度耦合的2对节点、不考虑刚度耦合的2对节点、单对节点的动力学模型，各模型预测结果如图4和表5所示。

对比各模型预测与实测结果发现：单对节点动力学模型在低频段容易丢失计算模态和出现多余计算模态，低频区的预测效果较差；含耦合关系的2对节点动力学有限元分析预测结果与实验实测结果的前11阶模态固有频率具有良好的一致性，且并未出现多余和丢失计算模态情况。其中，最小误差为0.21%，最大误差为3.84%，平均误差仅为1.6%；而不含耦合关系的2对节点动力学有限元分析预测结果与实验结果对比发现，虽未出现多余和丢失计算模态，但前11阶模态固有频率误差较大，平均误差高达17%。由此可知，为使有限元分析模型具有更好的预测结果，应考虑多对节点动力学建模方式及构建刚度信息间的耦合函数关系。

有限元分析模态振型与实验测得振型对比如图5所示。结果显示，两者振型吻合度较高，论证了方法的可行性及考虑刚度间耦合关系的必要性。

4 结　论

提出了基于响应面函数的螺栓结合部刚度与动力学性能指标—模态固有频率关系的建模方法，利用CCD实验设计方案，得到显现二者间的关联关系响应面函数模型，并结合遗传算法与二次规划寻优算法快速精准的实现螺栓结合部动刚度参数的辨识。

通过响应函数的选型，宏观上对比分析了结合部刚度间的耦合关系对结合部动力学参数辨识精度的影响。通过结合部单、双对节点的有限元动力学建模方式，对比分析了两种建模方式对结合部动力学性能的影响。由此可知节点对数的选择及刚度间的耦合关系对结合部动力学的研究有着重要价值与意义。

实例结果表明，为使响应面函数的辨识精度更高，在响应函数选型时应考虑含交叉关系的二次函数模型，充分考虑到刚度参数间的耦合及有限元建模方式对结合部动力学特性的影响。