基于NURBS啮合线的高效双螺杆转子自适应设计方法

Adaptive Design Method of Efficient Twin Screw Rotors Based on NURBS Mesh Lines

摘要

为解决压缩机转子型线设计过程繁琐以及性能难以评判的问题，对双螺杆压缩机转子型线的反向设计进行了研究。利用NURBS曲线对转子啮合线进行分段拟合，求解出各段啮合线的控制点，并建立转子型线的优化模型，利用单纯形法调整啮合线控制点坐标实现对转子型线的自适应优化设计。将优化前后的转子型线与原始型线比较分析并进行流体仿真，结果显示，优化后的转子型线综合性能提高了11.38%，优化后的压缩机转子表面压力数值在三种工况下分别提高了5.16%、5.00%、4.93%，证明了双螺杆压缩机转子型线自适应优化设计的有效性。

Abstract

To address the issues of cumbersome design processes and challenging performance evaluation in compressor rotor profile design， a reverse design approach for twin-screw compressor rotor profiles was investigated. The rotor meshing line was segmented and fitted using NURBS curves， the control points of each meshing line segment were solved， and an optimization model for the rotor profiles was established. An adaptive optimization design of the rotor profile was achieved by adjusting the coordinates of the meshing line control points using the simplex method. A comparative analysis between the optimized rotor profiles and the original profiles was conducted， along with fluid simulation. The results demonstrate that the comprehensive performance of the optimized rotor profiles is improved by 11.38%. Under three operating conditions， the surface pressure values of the optimized compressor rotors are increased by 5.16%， 5.00%， and 4.93%， respectively， which verify the effectiveness of the adaptive optimization design for twin-screw compressor rotor profiles.

Graphical abstract

关键词

转子型线 / 反向设计 / 啮合线 / 优化设计 / 非均匀有理B样条曲线

Key words

rotor profile / reverse design / meshing line / optimization design / NURBS curve

Author summay

耿宇航，男，2000年生，硕士研究生。研究方向为双螺杆压缩机转子型线的设计与优化。E-mail：[e-mail]

何雪明^*（通信作者），男，1966年生，教授。研究方向为自由曲面、曲线、CAD、CAM。E-mail： [e-mail]。

0 引言

双螺杆压缩机的核心部件是一对相互啮合的螺杆转子，转子型线的设计对压缩机效率、稳定性和耐用性等起到关键作用，所以深入研究和优化转子型线对提升压缩机的性能至关重要。转子型线设计方法有正向和反向两种。正向设计方法包括转子法和齿条法，其原理都是先定义阴阳转子中一条型线，再通过包络条件式和坐标转换关系推导出另一条型线，现有的大多数转子型线如SRM型线、日立型线等都是通过转子法进行设计的。徐健等^［1］基于齿轮啮合原理提出了齿廓法线法，推导了多种类型齿曲线的包络条件式，用矩阵方程描述坐标变换，提高了设计效率。卢家伦等^［2］利用三次参数样条拟合计算共轭型线，解决了传统二次曲线拼接出现奇异点的问题。YAKUPOV等^［3］将数值方法运用到转子型线的设计与分析中，使用Romberg方法计算接触线长度等参数，解决了传统解析法计算型线几何参数难以满足复杂型线分析需求的问题，提高了计算精度和效率。PATIL等^［4］运用路径同伦法为转子型线的设计提供了更科学的曲线选择，增强了曲线设计的灵活性。任焕梅等^［5］用三次B样条曲线表达转子型线，以接触线为优化目标对转子型线进行优化设计，并对优化前后转子型线进行流场仿真验证。虽然正向设计的研究较多，理论较为成熟，但正向设计无法直接评估型线性能，需设计完成后计算啮合线方程才能进行性能评估，若性能不符合要求修改型线会导致大量重复计算。

反向设计即啮合线法由ZAYTSEV等^［6-7］提出，通过定义啮合线设计转子型线，解决了正向设计难以预测性能的问题。徐健等^［8］推导了啮合线分别为点、直线、圆弧时对应的转子齿曲线方程。何雪明等^［9］详细分析了啮合线与转子型线之间关系，并给出了啮合线连续性的要求及求解过程。在啮合线定义完成后，其性能可能不满足预期要求，还需对其修改调整进行优化以达到性能要求。WU等^［10］提出使用三次样条曲线拟合啮合线，并利用二次方程调整啮合线以优化转子型线。HOANG等^［11］使用直线和椭圆曲线显式方程预定义啮合线生成不同齿数的转子型线并进行优化。

由于现有型线的啮合线大多为不规则曲线，而非均匀有理B样条（NURBS）曲线具有对不同曲线进行统一表达和局部修改的特点，且相比于其他自由曲线，NURBS曲线具有更高的精确性和更灵活的曲线修改能力，因此本文提出一种基于NURBS啮合线设计转子型线的方法。针对现有型线，利用NURBS曲线对其啮合线进行拟合，建立转子型线的优化模型，调整曲线控制点对啮合线进行修改，得到性能更优的转子型线，实现双螺杆压缩机转子型线的自适应寻优。

1 转子型线的反向设计

1.1 坐标变换关系

转子型线的反向设计需建立图1所示的5个坐标系。固接于机壳上的啮合线坐标系为O₀x₀y₀，阳阴转子的静坐标系分别为O₁X₁Y₁、O₂X₂Y₂，随阳阴转子转动的动坐标系分别为O₁x₁y₁、O₂x₂y₂。R₁、R₂为阳阴转子节圆半径，R_a1为阳转子齿顶圆半径，R_f2为阴转子齿根圆半径。H为阳阴转子中心距，满足关系式H=R₁+R₂。

双螺杆压缩机的阳阴转子可看作一对相互啮合的斜齿轮，转子之间是定传动比啮合关系，即传动比i满足

i=\frac{{\phi }_{\mathrm{2}}}{{\phi }_{\mathrm{1}}}=\frac{n_{\mathrm{2}}}{n_{\mathrm{1}}}=\frac{{\omega }_{\mathrm{2}}}{{\omega }_{\mathrm{1}}}=\frac{R_{\mathrm{1}}}{R_{\mathrm{2}}}=\frac{z_{\mathrm{1}}}{z_{\mathrm{2}}}（1）

式中：φ为阳阴转子旋转角度；n为转速；ω为角速度；R为节圆半径；z为齿数；下标1和2分别表示该参数为阳转子或阴转子。

令k=1+i，动坐标系O₁x₁y₁与O₀x₀y₀之间满足以下变换关系：

\left.\begin{array}{l}{x}_{\mathrm{1}}=x_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{1}}+y_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{1}}-R_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{1}}\\ y_{\mathrm{1}}=-x_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{1}}+y_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{1}}+R_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{1}}\end{array}\right\}（2）

动坐标系O₂x₂y₂与O₀x₀y₀之间满足

\left.\begin{array}{l}{x}_{\mathrm{2}}=x_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{2}}-y_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{2}}+R_{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{2}}\\ y_{\mathrm{2}}=x_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{2}}+y_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{2}}+R_{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{2}}\end{array}\right\}（3）

联立式（2）和式（3），可得坐标系O₂x₂y₂与O₁x₁y₁之间的变换关系：

\left.\begin{array}{l}{x}_{\mathrm{2}}=x_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}k{\phi }_{\mathrm{1}}-y_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}k{\phi }_{\mathrm{1}}+H\mathrm{c}\mathrm{o}\mathrm{s}i{\phi }_{\mathrm{1}}\\ y_{\mathrm{2}}=x_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}k{\phi }_{\mathrm{1}}+y\mathrm{c}\mathrm{o}\mathrm{s}k{\phi }_{\mathrm{1}}+H\mathrm{s}\mathrm{i}\mathrm{n}i{\phi }_{\mathrm{1}}\end{array}\right\}（4）

1.2 转子型线反向设计流程

转子型线反向设计是先定义啮合线的方程，再根据坐标转换关系和包络条件求解出相应的阴阳转子型线，反向设计流程见图2。

设某段啮合线的参数方程为

\begin{array}{c}{x}_{\mathrm{0}}=x_{\mathrm{0}}(t)\\ y_{\mathrm{0}}=y_{\mathrm{0}}(t)\end{array}t_{\mathrm{a}}\le t\le t_{\mathrm{b}}（5）

式中：t为啮合线方程参数；t_a、t_b为参数上下限。

将式（5）代入式（2）和式（3）可求得对应阳阴转子型线段的参数方程如下：

x_{\mathrm{1}}=x_{\mathrm{1}}(t,{\phi }_{\mathrm{1}})y_{\mathrm{1}}=y_{\mathrm{1}}(t,{\phi }_{\mathrm{1}})（6）

x_{\mathrm{2}}=x_{\mathrm{2}}(t,{\phi }_{\mathrm{2}})y_{\mathrm{2}}=y_{\mathrm{2}}(t,{\phi }_{\mathrm{2}})（7）

求解出阳阴转子的参数方程后还需求解包络条件式，即转角φ关于参数t的方程φ = φ（t）。根据啮合定律，双螺杆压缩机工作时阴阳转子在啮合点处的相对速度相等，相对速度方向和公切线方向相同，则相对速度方向和法线方向相互垂直，即

\mathit{n}·\mathit{v}=\mathrm{0}（8）

在阴转子坐标系下，啮合点处速度向量和法向量如下：

\mathit{v}=(\frac{\partial x_{\mathrm{2}}}{\partial {\phi }_{\mathrm{1}}},\frac{\partial y_{\mathrm{2}}}{\partial {\phi }_{\mathrm{1}}})\mathit{n}=(\frac{\partial y_{\mathrm{2}}}{\partial t},-\frac{\partial x_{\mathrm{2}}}{\partial t})（9）

将式（3）、式（9）代入式（8）可得φ₁与t的关系式：

k{y}_{\mathrm{0}}(t)\frac{\partial y_{\mathrm{0}}(t)}{\partial t}+k{x}_{\mathrm{0}}(t)\frac{\partial x_{\mathrm{0}}(t)}{\partial t}+k{R}_{\mathrm{1}}{y}_{\mathrm{0}}(t){\dot{\phi }}_{\mathrm{1}}=\mathrm{0}（10）

再对式（10）积分可求出包络条件式：

{\phi }_{\mathrm{1}}=-{\int }_{t_{\mathrm{a}}}^{t_{\mathrm{b}}}\frac{y_{\mathrm{0}}{\dot{y}}_{\mathrm{0}}^{}+x_{\mathrm{0}}{\dot{x}}_{\mathrm{0}}^{}}{R_{\mathrm{1}}{y}_{\mathrm{0}}}\mathrm{d}t+{\phi }_{\mathrm{0}}（11）

1.3 NURBS啮合线反向设计型线

由于啮合线种类繁多，难以用初等函数对其进行统一表达，可采用NURBS曲线来构造啮合线，且NURBS曲线具有局部修改的优势，可提高型线设计的效率和灵活性。

一条p次NURBS曲线可由4个参数进行定义：p次B样条基函数N_l，p （u）、权因子w_l 、节点矢量 U 及控制点P_l。采用有理分式的形式可表达为

C(u)=\frac{{\displaystyle \sum _{l=\mathrm{0}}^n}{N}_{l,p}(u)w_l{P}_l}{{\displaystyle \sum _{l=\mathrm{0}}^n}{N}_{l,p}(u)w_l}u_{\mathrm{m}\mathrm{i}\mathrm{n}}\le u\le u_{\mathrm{m}\mathrm{a}\mathrm{x}}（12）

采用逼近法构造一条NURBS曲线实现对啮合线的分段拟合。先在原始啮合线上选取一组描述曲线形状的数据点即型值点q_l （l=0，1，…，m），利用弦长参数化的方法对型值点进行参数化并归一化计算出参数值ū_l，求解节点矢量 U =（u₀，u₁，…，u_n+p+1），节点矢量两端的外节点满足p+1次重复度，即u₀=u₁=…=u_p =0，u_n+1=u_n+2=…=u_n+p+1=1。内节点通过下式计算：

u_{j+p}=(\mathrm{1}-t){\overline{u}}_{l-\mathrm{1}}+t{\overline{u}}_lc=\frac{m}{n-p}（13）

l=\mathrm{i}\mathrm{n}\mathrm{t}(jc)t=jc-lj=\mathrm{1,2},\cdot \cdot \cdot ,n-p

式中：m为数据点个数；n为控制点个数；p为曲线次数；{\overline{u}}_l为参数值，u为节点矢量中的节点；int（jc）表示小于等于jc的最大正整数。

再通过下式计算出基函数：

\left.\begin{array}{l}{N}_{l,\mathrm{0}}(u)=\left\{\begin{array}{l}\mathrm{1}{u}_l\le u\le u_{l+\mathrm{1}}\\ \mathrm{0}\mathrm{其他}\end{array}\right.p=\mathrm{0}\\ N_{l,p}(u)=\frac{u-u_l}{u_{l+p}-u_l}{N}_{l,p-\mathrm{1}}(u)+\\ \frac{u_{l+p+\mathrm{1}}-u}{u_{l+p+\mathrm{1}}-u_{l+\mathrm{1}}}{N}_{l+\mathrm{1},p-\mathrm{1}}(u)p>\mathrm{0}\end{array}\right\}（14）

为避免非线性问题，设定权因子w = 1，由最小二乘理论最小化构造的NURBS曲线与原始曲线数据之间的误差值^［12］可得求解各段曲线控制点的关系式：

\mathit{P}=({\mathit{N}}^{\mathrm{T}}{\mathit{N})}^{-\mathrm{1}}\mathit{R}（15）

式中： P 为控制点矩阵； N 为k次基函数矩阵； R 为误差值矩阵。

误差值r_l 和误差值矩阵 R 由下式计算：

\left.\begin{array}{l}{r}_l=q_l-N_{\mathrm{0},p}(w_l)q_{\mathrm{0}}-N_{n,p}(w_l)q_{m}l=\mathrm{1,2},\dots ,m-\mathrm{1}\\ R_l={\displaystyle \sum _{j=\mathrm{0}}^{m-\mathrm{1}}}(N_{l,p}{w}_j{r}_j)l=\mathrm{0,1},\dots ,n-\mathrm{1}\end{array}\right\}（16）

拟合后的NURBS啮合线可由下式表示：

\begin{array}{c}{x}_{\mathrm{0}}=C_x(u)\\ y_{\mathrm{0}}=C_y(u)\end{array}\mathrm{0}\le u\le \mathrm{1}（17）

\left.\begin{array}{c}\\ \end{array}\right\}

将式（17）代入包络条件式（式（11））可得NURBS啮合线的包络条件式：

{\phi }_{\mathrm{1}}=-\frac{\mathrm{1}}{R_{\mathrm{1}}}{\int }_{\mathrm{0}}^u\frac{C_y(u){\dot{C}}_y^{}(u)+C_x(u){\dot{C}}_x^{}(u)}{C_y(u)}\mathrm{d}t+{\phi }_{\mathrm{0}}（18）

φ₀为每段啮合线啮合起始角度，啮合线第一段曲线取φ₀=0。求解φ₁还需对式（18）进行积分，由于被积函数无法用初等函数表达，且为了利于编程计算，可用数值积分方法中龙贝格积分公式对其进行积分计算^［13］，再将NURBS啮合线和包络条件式代入式（6）、式（7），即可求出对应的阴阳转子型线方程。

2 转子型线的优化

基于NURBS理论对转子型线进行反向设计和自适应优化。先选取一条原始啮合线，对各段啮合线提取型值点，通过逼近的方法对啮合线进行拟合，构造出一条NURBS啮合线，建立转子型线的优化模型，调整NURBS啮合线的控制点，对NURBS啮合线进行修改，经过多次迭代，设计出一条性能最优的转子型线。

2.1 原始啮合线的计算与分段

本文以第三代不对称转子型线中的复盛型线为例进行反向设计及优化，如图3所示。复盛型线由四段齿曲线组成，如表1所示。其中，阴阳转子齿数比为6∶5，中心距为110 mm。阳转子齿顶圆半径R_a1=77 mm，阴转子齿顶圆半径R_a2=61 mm。

由表1可以看出，复盛型线齿曲线由二次曲线圆弧和椭圆弧及其包络线组成。传统正向设计方法先定义二次曲线的方程，根据包络条件式求解出对应二次曲线的包络线和啮合线公式完成一条型线的设计，该条复盛型线的啮合线参数方程如下：

AB段齿曲线的啮合线参数方程

\left.\begin{array}{l}{x}_{\mathrm{0}}=\mathrm{59.4}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{a}\mathrm{b}}+\mathrm{1.6}\mathrm{c}\mathrm{o}\mathrm{s}(t-{\phi }_{\mathrm{a}\mathrm{b}})\\ y_{\mathrm{0}}=-\mathrm{59.4}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{a}\mathrm{b}}+\mathrm{1.6}\mathrm{s}\mathrm{i}\mathrm{n}(t-{\phi }_{\mathrm{a}\mathrm{b}})\\ {\phi }_{\mathrm{a}\mathrm{b}}(t)=t-\frac{\mathrm{\pi }}{\mathrm{180}}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{0.97}\mathrm{s}\mathrm{i}\mathrm{n}t)\\ \mathrm{0}\le t\le \mathrm{1.92}\mathrm{r}\mathrm{a}\mathrm{d}\end{array}\right\}（19）

BC段齿曲线的啮合线参数方程

\left.\begin{array}{l}{x}_{\mathrm{0}}=-\mathrm{67.8}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{b}\mathrm{c}}-\mathrm{9.2}\mathrm{c}\mathrm{o}\mathrm{s}(t-{\phi }_{\mathrm{b}\mathrm{c}})\\ y_{\mathrm{0}}=-\mathrm{67.8}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{b}\mathrm{c}}-\mathrm{9.2}\mathrm{s}\mathrm{i}\mathrm{n}(t+{\phi }_{\mathrm{b}\mathrm{c}})\\ {\phi }_{\mathrm{b}\mathrm{c}}(t)=\frac{\mathrm{\pi }}{\mathrm{180}}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{1.356}\mathrm{s}\mathrm{i}\mathrm{n}t)-t\\ \mathrm{0}\le t\le \mathrm{0.76}\mathrm{r}\mathrm{a}\mathrm{d}\end{array}\right\}（20）

CD段齿曲线的啮合线参数方程

\left.\begin{array}{l}{x}_{\mathrm{0}}=-\mathrm{35.26}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{c}\mathrm{d}}-\mathrm{41.74}\mathrm{c}\mathrm{o}\mathrm{s}t\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{c}\mathrm{d}}-\\ \mathrm{35.4648}\mathrm{s}\mathrm{i}\mathrm{n}t\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{c}\mathrm{d}}\\ y_{\mathrm{0}}=\mathrm{35.26}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{c}\mathrm{d}}+\mathrm{41.74}\mathrm{c}\mathrm{o}\mathrm{s}t\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{c}\mathrm{d}}+\\ \mathrm{35.4648}\mathrm{s}\mathrm{i}\mathrm{n}t\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{c}\mathrm{d}}\\ {\phi }_{\mathrm{c}\mathrm{d}}(t)=[\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}(\frac{\mathrm{0.24}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{2}t-\mathrm{1.68}\mathrm{s}\mathrm{i}\mathrm{n}t}{\sqrt[]{\mathrm{4.33}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}t+\mathrm{3.15}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}t}})-\\ \mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}(\mathrm{1.18}\mathrm{t}\mathrm{a}\mathrm{n}t)]\frac{\mathrm{\pi }}{\mathrm{180}}\\ \mathrm{0}\le t\le \mathrm{1.57}\mathrm{r}\mathrm{a}\mathrm{d}\end{array}\right\}（21）

DE段齿曲线的啮合线参数方程

\left.\begin{array}{l}{x}_{\mathrm{0}}=-\mathrm{52.636}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{d}\mathrm{e}}+\mathrm{3.636}\mathrm{c}\mathrm{o}\mathrm{s}(t-{\phi }_{\mathrm{d}\mathrm{e}})\\ y_{\mathrm{0}}=-\mathrm{52.636}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{d}\mathrm{e}}+\mathrm{3.636}\mathrm{s}\mathrm{i}\mathrm{n}({\phi }_{\mathrm{d}\mathrm{e}}-t)\\ {\phi }_{\mathrm{d}\mathrm{e}}(t)=t-\frac{\mathrm{\pi }}{\mathrm{180}}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{1.05}\mathrm{s}\mathrm{i}\mathrm{n}t)\\ \mathrm{0}\le t\le \mathrm{0.73}\mathrm{r}\mathrm{a}\mathrm{d}\end{array}\right\}（22）

采用反向设计方法直接对啮合线进行修改，使型线性能符合设计要求。为便于对啮合线进行NURBS曲线拟合以及自适应优化，将原始啮合线在特定极点处分为八段，如图4所示。以阴阳转子节圆切点B点为原点建立坐标系，将啮合线分为四个象限，其中一、三象限啮合线对应高压侧齿廓，二、四象限啮合线对应低压侧齿廓。从最右侧啮合线与x轴交点A点处开始，沿图中箭头方向分为A-M-B-C-D-E-B-N-A八段。A点为阳转子齿根圆与阴转子齿顶圆的切点，D点为阳转子齿顶圆与阴转子齿根圆的切点，改变A、D点坐标位置可改变阴阳转子齿顶圆与齿根圆的直径。C点为啮合线上的尖点，是转子型线上一个啮合周期的开始。M、N、E点分别为AB、BA、DB各段曲线的最高点或最低点。

2.2 转子型线的优化模型

优异的螺杆转子应具有较大的齿间容积，且齿间容积之间应具有良好的密封性。齿间容积是压缩机进行气体压缩的工作腔，影响压缩机的工作效率和性能，齿间容积主要受转子型线上三个几何参数的影响：接触线长度L、泄漏三角形面积S和面积利用系数M。根据转子型线的设计准则，转子型线应保证接触线长度较短、泄漏三角形面积较小、面积利用系数较大^［14］。由于三个参数之间相互影响，即接触线长度受到转子型线齿曲线形状的影响，齿曲线形状的改变又会对面积利用系数M的值产生影响，而面积利用系数与泄漏三角形面积之间存在相互制约的关系，难以只利用其中一个参数作为型线的评判指标，所以需建立一个综合评价指标来评估三个参数对转子型线性能的影响^［15］。引用推出的压缩机转子型线综合性能指标K作为型线优化的评判指标^［16-17］，K通过多目标决策中的模糊综合评价法和层次分析法分析各参数之间的关系和在综合指标中所占权重确定，K值越小，型线的综合性能越好，K的表达式如下：

K=\frac{\mathrm{1}}{\mathrm{60}}(\mathrm{0.1741}\sqrt[]{L}+\frac{\mathrm{0.7225}S}{\mathrm{0.1033}M})（23）

设计啮合线时，还需满足三个限制条件。

1）几何约束。啮合线必须位于阴阳转子型线齿顶圆的重叠区域，否则转子之间会发生干涉。

2）齿曲线干涉约束。转子型线单个齿曲线的圆心角θ\le\frac{\mathrm{2}\mathrm{\pi }}{z}。啮合线设置不合理导致θ>\frac{\mathrm{2}\mathrm{\pi }}{z}时，转子型线的齿曲线的末端与下一个齿的齿曲线首端会出现干涉。

3）底切限制。啮合线设置不合理会使得啮合线上点的切线方向和该点与阴阳转子轴心连线方向垂直，导致型线上出现奇异点。由于

\left.\begin{array}{l}\mathit{\tau }·{\mathit{n}}_{\mathrm{1}}=\mathrm{0}\mathit{\tau }·{\mathit{n}}_{\mathrm{2}}=\mathrm{0}\\ \mathit{\tau }=({\dot{C}}_x^{}(u),{\dot{C}}_y^{}(u))\\ {\mathit{n}}_{\mathrm{1}}=(C_x(u)-R_{\mathrm{1}},C_y(u))\\ {\mathit{n}}_{\mathrm{2}}=(C_x(u)+R_{\mathrm{2}},C_y(u))\end{array}\right\}（24）

其中， τ 为啮合线上点的切线向量， n₁、 n₂为啮合线上点指向阴阳转子轴心的方向向量，那么为避免干涉产生，需满足参数值u≤u_uc^［10］，其中u_uc为发生底切的参数值。以各段NURBS啮合线的控制点坐标P_xl 、P_yl 为设计变量，以型线综合性能指标K值减小为优化目标，建立优化模型如下：

\mathrm{m}\mathrm{i}\mathrm{n}f=K=f(P_{xl},P_{yl})

\mathrm{s}.\mathrm{t}.g_{\mathrm{1}}=\sqrt[]{C_x(u)P_{xl}^{\mathrm{2}}+C_y(u)P_{yl}^{\mathrm{2}}}-R_{\mathrm{a}\mathrm{1}}\le \mathrm{0}

g_{\mathrm{2}}=\sqrt[]{C_x(u)P_{xl}^{\mathrm{2}}+C_y(u)P_{yl}^{\mathrm{2}}}-R_{\mathrm{a}\mathrm{2}}\le \mathrm{0}

g_{\mathrm{3}}=\theta (P_{xl},P_{yl})-\frac{\mathrm{2}\mathrm{\pi }}{z}\le \mathrm{0}

g_{\mathrm{4}}=u(P_{xl},P_{yl})-u_{\mathrm{u}\mathrm{c}}\le \mathrm{0}

式中：f为目标函数；g₁、g₂为啮合线的几何约束条件；g₃为齿曲线干涉约束条件；g₄为底切限制约束条件。

2.3 转子型线的自适应优化

对原始复盛型线的八段啮合线每段选取45个型值点，利用5次NURBS曲线对各段啮合线进行拟合，求解出每段啮合线的6个控制点。根据2.2节中建立的优化模型，通过单纯形法自适应调整NURBS啮合线控制点坐标，约束条件通过惩罚函数法纳入目标函数中，如下式所示：

F(P_{xl},P_{yl})=f(P_{xl},P_{yl})+r{\displaystyle \sum _{l=\mathrm{1}}^{\mathrm{5}}}\mathrm{m}\mathrm{a}\mathrm{x}(\mathrm{0},g_l(P_{xl},P_{yl}{))}^{\mathrm{2}}（25）

其中，r为惩罚因子。当约束条件满足时，式（25）右侧障碍项的值为0；违反约束条件时，障碍项值为正值，使目标函数值增大从而停止优化迭代。 单纯形法属于无约束优化方法中的一种，通过目标函数值的信息来构建搜索方向。单纯形即n维空间中的n+1个顶点的多面体，对啮合线控制点坐标的调整是在二维平面上进行的，此时建立的单纯形即为平面上的三角形，先计算目标函数在若干点处函数值构建初始单纯形，再通过不断的迭代构建新的单纯形，以搜索目标函数的极小点，其具体过程如下：

1）建立初始单纯形，计算初始点x₀处的函数值f₀，根据设定步长，选取不在同一直线上的另外两个顶点x₁、x₂并计算函数值f₁、f₂。

2）比较三个顶点函数值大小确定最好点x_L，最差点x_H和次差点x_G。其中，令f_L=f（x_L）、f_H=f（x_H）、f_G=f（x_G）。

3）检验是否满足收敛准则|（f_L-f_H）/f_L|ɛ，ɛ为收敛精度值。若满足，则f_min=f_L，结束搜索。否则进行下一步。

4）计算x_L和x_G两点的重心x₃，x₃=0.5（x₀+x₁+x₂-x_H），反射点x₄=2x₃-x_H及反射点的函数值f₄=f（x₄）。①当f_L≤f₄f_G时，进行反射操作，用x₄代替x_H，f₄代替f_H，构成新的单纯形，然后回到第二步进行循环。②当f₄f_L时，进行扩张操作，求取扩张点x₅=x₃+α（x₄-x₃）及其函数值f₅=f（x₅）。其中，α取1.2~2.0，使f₅f₄，用x₅代替x_H，f₅代替f_H，构成新的单纯形，返回步骤2）进行循环。③当f₄≥f_G时进行收缩操作。求收缩点x₆=x₃+β（x₄-x₃）和收缩点函数值f₆=f（x₆），其中，β一般取0.5~1.0，若f₆f_H，用x₆代替x_H得到新的单纯形，返回到步骤2）进行循环。若f₆≥f_H，以最好点x_L为中心缩边，使单纯形的其他两个顶点向x_L移动一半距离，得到新的单纯形，然后返回到步骤2）进行循环。

每段啮合线的首尾控制点坐标保持不变，对其他控制点（P_xl，P_yl ）通过以上步骤在0.010、0.005和0.001的步长下迭代进行自适应调整，计算出性能最优的一条啮合线。优化前后三个参数值和型线综合性能指标K的变化如表2所示，优化前后啮合线变化如图5所示，型线齿廓的变化如图6所示。与原始型线相比，优化后的转子型线接触线长度减小了0.3725 mm，泄漏三角形面积减小了1.274 mm²，面积利用系数增大了0.0064，型线综合性能指标K减小了0.388，型线综合性能提高了11.38%。

3 流体动力学仿真与实验验证

3.1 流体域模型与仿真参数设置

对优化前后的转子在UG三维软件中进行建模，并建立相应的流体域模型（图7）。通过ANSA网格划分软件对阴阳转子进行面网格划分，导入STAR CCM+仿真软件进行流体域面网格和体网格划分，如图8所示，进气腔和排气腔网格类型设置为多面体网格，压缩腔网格类型设置为切割体网格，以减少模型网格数量，提高计算效率。双螺杆压缩机的CFD物理模型设置为以下选项：介质设为理想气体，时间模型选择隐式非稳态模型，湍流模型选择Realizable k-ε模型，流动模型选择耦合流动模型。

在阳转子转速分别为2940 r/min、3480 r/min和4020 r/min，压缩机进气压力均为大气压力0.101 MPa，排气端压力均为0.9 MPa的三组工况下对优化前后的转子进行流体动力学仿真。其中，第一组工况的仿真参数设置与原始复盛转子压缩机的实验工况一致，第二、第三组工况的阳转子转速依次提高540 r/min。

3.2 仿真结果对比分析

在STAR CCM+仿真软件中对优化前后的转子在三组工况下进行流体动力学仿真计算，计算结果如图9~图11所示。由仿真云图可以看出，转子表面压力均在排气端达到最大值，在三个工况下优化后的转子表面压力数值均得到了不同程度的增大。对转子型线优化后，阳转子转速n为2940 r/min时，转子表面最大压力p_m从0.920 03 MPa提高到0.967 48 MPa，提高了5.16%；阳转子转速n为3480 r/min时，转子表面最大压力p_m从0.917 63 MPa增大到0.963 47 MPa，提高了5.00%；阳转子转速n为4020 r/min时，转子表面最大压力p_m从0.913 21 MPa增大到0.958 22 MPa，提高了4.93%。三种工况下转子表面压力值的增大速率也得到了提高，说明转子齿间容积内气体的泄漏有所减少，证明在对转子型线进行自适应寻优后，转子齿间容积的密封性能得到了提高。

由图12可以看出，随着转子转速提高，转子表面压力值逐渐下降。这是因为压缩机转速越高，内部流量也越大，压缩机内部能量的损失也会提高，造成内部气体的压力值随之减小。

3.3 实验验证

为验证仿真结果的有效性和可靠性，搭建原始复盛型线的双螺杆压缩机实验平台。实验设备包含双螺杆压缩机整机、数据采集卡、压力传感器和计算机等，如图13所示。在压缩机机壳上钻取4个直径为6 mm、深度为10 mm、间距为40 mm的螺纹孔安装传感器， 在压缩机阳转子转速2940 r/min、排气压力0.9 MPa的工况下，采集压缩机工作过程中的压力值，监测点布置示意见图14。

根据压缩机上设置的监测点位置，对复盛型线优化前后的转子进行流体仿真时流体域的相应位置设置探针，如图15所示，记录4个监测点的压力变化。

压缩机采集的数据和仿真记录数据对比如图16所示，可以看出，监测点1~4处压力的实验值和转子型线优化前后的仿真值均呈周期性波动，实验测量值与原始型线的仿真值重合较好，且优化后的转子型线在各监测点的压力值与优化前的相比得到了一定的提高，验证了流体动力学仿真的有效性和结果的可靠性。

4 结论

本文提出一种基于NURBS啮合线设计转子型线的方法，以复盛型线为例进行优化设计，建立优化模型，调整NURBS啮合线控制点对啮合线进行修改，设计出一条性能优异的转子型线。对优化前后转子型线进行对比并进行流体仿真验证，得到以下结论：

1）对自适应优化前后转子型线的几何参数进行对比，优化后的型线接触长度值减小了0.3725 mm、泄漏三角形面积值减小了1.274 mm²，面积利用系数值增大了0.0046，型线的综合性能提高了11.38%。

2）通过STAR CCM+仿真软件在三种工况下对优化前后的转子型线进行流体仿真验证，优化后转子表面压力分别增大了5.16%、5.00%、4.93%，证明优化后的压缩机转子齿间容积密封性更好，验证了转子型线自适应优化的有效性。最后搭建了实验平台，将实验结果和仿真结果进行对比，验证了流体动力学仿真的可靠性和准确性。

参考文献

[1] 徐健， 余宾宴， 余小玲， 等. 螺杆压缩机转子型线设计方法［J］. 压缩机技术， 2012 （2）： 1-6.

[2] XU Jian， YU Binyan， YU Xiaoling， et al. Design Method for Rotor Profile of Screw Compressor［J］. Compressor Technology， 2012 （2）： 1-6

[3] 卢家伦， 田雅芬， 邢子文. 基于离散数据的双螺杆压缩机共轭转子型线计算［J］.压缩机技术， 2019 （3）： 1-5.

[4] LU Jialun， TIAN Yafen， XING Ziwen. Conjugate Rotor Profile Calculation of Twin-Screw Compressor Based on Discrete Data［J］. Compressor Technology， 2019 （3）： 1-5.

[5] YAKUPOV R R， MUSTAFIN T N， KHAMIDULLIN M S， et al. Numerical Methods for Calculation of Geometrical Criteria for the Assessment of Screw Compressor Rotor Profile［J］. AIP Conference Proceedings， 2019， 2141（1）： 030059.

[6] PATIL S， STOSIC N， KOVACEVIC A， et al. Application of Path Homotopy in Twin Screw Compressor Rotor Profile Design［C］∥IOP Conference Series： Materials Science and Engineering.Chicogo 2021， 1180： 012004.

[7] 任焕梅， 葛荣雨， 陈英雷. 双螺杆压缩机转子型线设计与优化［J］. 机械设计， 2021， 38（4）： 102-106.

[8] REN Huanmei， GE Rongyu， CHEN Yinglei. Design and Optimization of Rotor Profile for Twin-screw Compressor［J］. Mechanical Design， 2021， 38（4）： 102-106.

[9] ZAYTSEV D， INFANTE FERREIRA C A. Profile Generation Method for Twin Screw Compressor Rotors Based on the Meshing Line［J］. International Journal of Refrigeration， 2005， 28（5）： 744-755.

[10] ZAYTSEV D， INFANTE FERREIRA C A. Screw Compressor Rotors Profile Generation Method Based on Pre-defined Meshing Line［C］∥3rd International Compressor Technique Conference.Xi'an. 2001：117-127.

[11] 徐健，周恩民，许岭松，等.基于啮合线法的螺杆转子型线设计方法［J］. 压缩机技术， 2013 （1）： 1-7.

[12] XU Jian， ZHOU Enmin， XU Lingsong， et al. Design Method of Screw Rotor Profile Based on Meshing Line Method［J］. Compressor Technology， 2013 （1）： 1-7

[13] 何雪明， 姜振钢， 武美萍， 等. 双螺杆压缩机转子啮合线与型线关系研究［J］. 中国机械工程， 2017， 28（7）： 763-770.

[14] HE Xueming， JIANG Zhengang， WU Meiping， et al. Research on the Relationship Between the Meshing Line and Profile of the Rotor of Twin-screw Compressor［J］. China Mechanical Engineering， 2017， 28（7）： 763-770.

[15] WU Y R， FONG Z H. Improved Rotor Profiling Based on the Arbitrary Sealing Line for Twin-screw Compressors［J］. Mechanism and Machine Theory： Dynamics of Machine Systems Gears and Power Trandmissions Robots and Manipulator Systems Computer-Aided Design Methods， 2008， 43（6）： 695-711.

[16] HOANG M T， TRAN T V. Methodology for Generating a General Profile of Asymmetrical Rotors of Twin-screw Vacuum Pumps with Point Meshing Feature by a Pre-defined Sealing Line［J］. Proceedings of the Institution of Mechanical Engineers， Part B. Journal of Engineering Manufacture， 2023， 237（10）： 1582-1595.

[17] PEIGL L. 非均匀有理B样条 ［M］. 2版.赵罡，译.北京： 清华大学出版社， 2010.

[18] PEIGL L. Non Uniform Rational B-splines ［M］. 2nd ed. ZHAO Gang，Trans. Beijing： Tsinghua University Press， 2010.

[19] 黄海楠， 何雪明， 张荣， 等. 基于NURBS曲线的螺杆转子型线反向设计研究［J］. 食品与机械， 2018， 34（3）： 93-97.

[20] HUANG Hainan， HE Xueming， ZHANG Rong， et al Research on Reverse Design of Screw Rotor Profile Based on NURBS Curve［J］. Food and Machinery， 2018， 34（3）： 93-97.

[21] 邢子文. 螺杆压缩机——理论、设计及应用［M］. 北京： 机械工业出版社， 2000.

[22] XING Ziwen. Screw Compressor —Theory， Design and Application［M］. Beijing： Machinery Industry Press， 2000.

[23] 何雪明， 胡蓉， 于嘉川， 等. 一种双螺杆转子型线综合性能判断及优化设计方法：CN202111412125.1 ［P］. 2023-08-15.

[24] HE Xueming， HU Rong， YU Jiachuan， et al. A Comprehensive Performance Evaluation and Optimization Design Method for Twin-screw Rotor Profile：CN202111412125.1 ［P］. 2023-08-15.

[25] 宋健忠. 基于NURBS曲线分段自动寻优的双螺杆压缩机转子啮合线的研究［D］.无锡： 江南大学， 2022.

[26] SONG Jianzhong. Research on Rotor Meshing Line of Twin-Screw Compressor Based on NURBS Curve Segmentation Automatic Optimization［D］. Wuxi： Jiangnan University， 2022.

[27] 宋健忠， 高宗， 何雪明， 等.基于综合性能与控制点的转子啮合线分段寻优［J］. 液压气动与密封， 2023， 43（7）： 12-16.

[28] SONG Jianzhong， GAO Zong， HE Xueming， et al. Segmented Optimization of Rotor Meshing Line Based on Comprehensive Performance and Control Points［J］. Hydraulic Pneumatic and Sealing， 2023， 43 （7）： 12-16.

基金资助

国家自然科学基金(51275210)

国家自然科学基金(51975251)

江苏省食品先进制造装备技术重点实验室自主研究课题(FMZ2018Y2)