0 引言双螺杆压缩机具有可靠性高、适应性强、操作和维护方便以及动态平衡良好等多种优点，已经广泛应用于冶金、食品、化工等多个行业。双螺杆压缩机的不断发展离不开国内外学者持之以恒的研究。张丽梅等[1]针对型线数据为离散点的情况，给出了一种计算双螺杆压缩机齿间面积的有效方法。严迪[2]利用流体动力学、动网格技术、啮合理论、微分几何学、多相流理论、空泡动力学理论、数值方法等，深入分析了泵腔三维瞬态流场特性，研究探讨了其内部工作原理及其性能影响因素。郭高锦[3]以双螺杆压缩机转子反向设计及动力学仿真作为切入点，在此基础上开展了深入的研究。Rane等[4]研制了关键帧重新网格化技术，以不同的时间步长给计算流体动力学（Computational Fluid Dynamics，CFD）求解器预生成非结构化网格，采用多种网格变形方法对可逆定热压缩过程成功进行了模拟。Mustafin等[5-6]提出了可利用数值方法计算工况下转子温度场的仿真模型，分析了螺杆压缩机结构参数及工况对螺杆转子温度场的影响。Buckney等[7]给出了一个在原始和修正间隙下升高排气温度来测试压缩机的实例研究，预测可能的热变形和间隙产生的影响，对原始和修正的间隙设计进行评估，并讨论了由于修改而产生的性能损失。然而，典型型线在综合性能的表现方面存在着一定的不足，且目前也没有一个较为明确的方法来对典型型线进行优化。另外，在典型型线调整后对压缩机综合性能的影响方面，目前的研究还较为欠缺，相应的调整缺乏一定的理论依据。本文从转子型线的正向设计出发，结合一种基于综合性能自动寻优的方法，得到优化后的第三代典型转子型线，并对优化后的转子型线进行了仿真。仿真结果表明，自动寻优后的转子型线具有较好的密封性能。1 转子型线的正向设计研究1.1　坐标转换设计螺杆转子型线必须要建立相应的坐标系，如图1所示。其中，O1X1Y1和O2X2Y2均为转子的静坐标系，前者为阳转子静坐标系，后者为阴转子静坐标系；相应的动坐标系分别为O1x1y1和O2x2y2。其中，阳转子绕O1点逆时针旋转，阴转子绕O2点顺时针旋转。10.16578/j.issn.1004.2539.2023.11.008.F001图1转子法中坐标系及其关系Fig. 1Coordinate system and its relation in rotor method坐标系统中一共有4种坐标系，它们之间的转换关系为X1=x1cosφ1+y1sinφ1Y1=-x1sinφ1+y1cosφ1 （1）X2=x2cos(iφ1)+y2sin(iφ1)Y2=-x2sin(iφ1)+y2cos(iφ1) （2）X1=H-X2Y1=Y2 （3）式中，i=φ2/φ1；H为中心距。整理式（1）~式（3）可得x1=-x2cos(kφ1)-y2sin(kφ1)+Hcosφ1y1=-x2sin(kφ1)+y2cos(kφ1)+Hsinφ1 （4）式中，k=1+i。1.2　型线以及啮合线求解螺杆转子的型线通常由多段齿曲线首尾相连组成，在设计时，往往是先确定某转子的一段转子型线后求其共轭曲线。当已知某转子上的某组成齿曲线段方程，求解另一个转子上与之共轭的曲线方程时，可运用解析包络法和齿廓法线法。现对齿廓法线法做详细介绍。在某一曲线坐标系中，阴、阳转子接触点的公法线通过瞬时回转中心，根据斜率方程可得Y1-y1X1-x1=Ny1Nx1 （5）接触点的坐标值为（x1，y1），瞬时回转中心的坐标值为（X1，Y1），二者都位于阳转子的动坐标系内。将这两点相连，连线的矢量方向和接触点的法向方向一致。Nx1和Ny1分别为接触点法向量的水平分量和垂直分量。易求得在阳转子动坐标系中，瞬时回转中心的坐标为(-rj1cos φ1，rjlsin φ1）。其中，rj1为阳转子的节圆半径。而接触点的坐标值可以通过图2所示某段动坐标系中的椭圆齿曲线说明[8]，有10.16578/j.issn.1004.2539.2023.11.008.F002图2齿廓法线法在椭圆中的应用Fig. 2Application of tooth profile normal method in ellipser1(t)=-(rt1-a)-acostbsint1 （6）式中，rt1为齿顶圆半径；a为椭圆的长半轴长；b为椭圆的短半轴长。其中，接触点的坐标值为[-(rt1-a)-acost，bsint]。将接触点坐标值、瞬时回转中心坐标值以及接触点法向量的水平分量和垂直分量代入式（5），整理可得阳转子动坐标系旋转角度和参数t之间的关系为φ1(t)=arcsinb2-a22sin(2t)-a(rtl-a)sintrj1(asint)2+(bcost)2+arctan(abtant) （7）式（7）为该段椭圆齿曲线的啮合线方程，其对应的椭圆包络线表达式为r2(t)=[-(rt1-a)-acost]cos[(1+i)φ1(t)]-bsintsin[(1+i)φ1(t)]+Hcos[φ1(t)i][-(rt1-a)-acost]sin[(1+i)φ1(t)]+bsintcos[(1+i)φ1(t)]+Hsin[φ1(t)i]1 （8）将长半轴a和短半轴b替换成圆弧半径r，即可得到圆弧齿曲线、圆弧包络线和对应的啮合线表达式，分别如式（9）~式（11）所示。r1(t)=-(rt1-r)-rcostrsint1 （9）r2(t)=-(rt1-r)cos[(1+i)φ1(t)]-rcos[t-(1+i)φ1(t)]+Hcos[φ1(t)i]-(rt1-r)sin[(1+i)φ1(t)]+rsin[t-(1+i)φ1(t)]+Hsin[φ1(t)i]1 （10）φ1(t)=t-arcsin(rt1-r)sintrj1 （11）2 基于NURBS的型线控制点求解2.1　NURBS曲线的定义NURBS（非均匀有理B样条）可分为有理基函数和有理分式等几种不同的表达形式，虽然侧重不同，但都是等价的。本文以有理分式为例进行说明。一条p次NURBS曲线的定义为C(u)=∑i=0nNi,p(u)wiPi∑i=0nNi,p(u)wi,a≤u≤b （12）式中，{Pi}为控制点，各控制点连接起来构成的多边形称为控制多边形；{wi}为权因子，{wi}的大小决定了曲线偏离控制多边形的程度，其中，首末权因子w0、wn0，其余wi≥0；Ni，p(u)为定义在节点矢量U上的p次B样条基函数，U中节点的分布非均匀且非周期，U={a，…，a︸p+1，up+1，…，um-p-1，b，…，b︸p+1}，一般情况下端节点分别取a=0，b=1。2.2　基函数的计算B样条基函数的定义方法有很多种。为计算方便，且易于编程计算实现，本文选用递推公式定义，分别如式（13）~式（14）所示。Ni,0=1,ui≤u≤ui+10,其他 （13）Ni,p(u)=u-uiui+p-uiNi,p-1(u)+ui+p+1-uui+p+1-ui+1Ni+1,p-1(u) （14）由式（13）、式（14）可知，0次基函数Ni，0(u)是一个阶梯函数，只在区间[ui，ui+1]上为1，其余为0；基函数大于0次时，基函数Ni，p(u)是两个低一次的基函数Ni，p-1(u)和Ni+1，p-1(u)的线性组合，由此递推关系可以得到任意次数的基函数。需要说明的是，在计算基函数之前要先指定次数p和节点矢量U。2.3　NURBS曲线的导矢无论是正向设计还是反向设计，都需要用到包络条件式，需要对曲线的参数方程进行求导；此外，在研究曲线间连续性问题时，也需要利用曲线的导矢进行分析。上述问题都涉及曲线的求导，因此，想用NURBS来表达转子型线和啮合线，就必须得到NURBS的求导方法。2.3.1 1阶导数假设A(u)=∑i=0nNi,p(u)wiPiK(u)=∑i=0nNi,p(u)wi （15）则式（12）可变为C(u)=A(u)K(u) （16）因此，根据分式求导公式对NURBS曲线C(u)求1阶导可得C'(u)=K(u)A'(u)-K'(u)A(u)K(u)2=K(u)A'(u)-K'(u)K(u)C(u)K(u)2=A'(u)-K'(u)C(u)K(u) （17）A'(u)和K'(u)可通过对基函数求导得到，有A'(u)=∑i=0nNi,p'(u)wiPiK'(u)=∑i=0nNi,p'(u)wi （18）而对于基函数的求导，可由下述递推公式求得Ni,p'(u)=pui+p-uiNi,p-1(u)+pui+p+1-ui+1Ni+1,p-1(u) （19）将式（17）~式（19）联立，即可求得NURBS曲线的1阶导矢。2.3.2 高阶导数若想得到NURBS曲线C（u）的更高阶导数，可以通过对A(u)利用牛顿-莱布尼茨公式进行求导，即A(k)(u)=[K(u)C(u)](k)=∑i=0kkiK(i)(u)C(k-i)(u)=K(u)C(k)(u)+∑i=1kkiK(i)(u)C(k-i)(u) （20）下面以2阶导数为例进行说明。由式（20）可得NURBS曲线的2阶导数为C(2)(u)=A(2)(u)-∑i=122iK(i)(u)C(2-i)(u)K(u)=A(2)(u)-2K(1)(u)C(1)(u)-K(2)(u)C(u)K(u) （21）A(u)与K(u)的2阶导数都可转换成关于基函数的2阶导数，即A"(u)=∑i=0nNi,p″(u)wiPiK″(u)=∑i=0nNi,p″(u)wi （22）式（22）中基函数的2阶导数可以由下式求得：Ni,p″(u)=pp-2[u-uiui+p-uiNi,p-1″(u)-ui+p+1-uui+p+1-ui+1Ni+1,p-1″(u)] （23）将式（21）~式（23）联立，即可求得C(u)的2阶导矢。2.4　型值点法求解控制点型值点法求解控制点，顾名思义，就是通过型线上的型值点求出转子型线上的控制点。要求出插值于qi的NURBS曲线，就必须确定节点矢量、控制点和权因子。对于权因子，默认取1，而节点矢量和控制点在下面进行详细说明。2.4.1 型值点的参数化首先，需要了解型值点的参数化。对于给定插值点qi，确定其对应的参数值ui的最常用方法是弦长参数化，其中，可以假定参数u∈[0，1]，即相邻参数ui的差值比与对应型值点之间的距离比相等。假定型值点qi的总弦长用d表示，由d=∑i=1nqi-qi-1 （24）易得u0=0,un=1,ui=ui-1+qi-qi-1d，i=1，2，…，n-1（25）2.4.2 节点矢量在确定了型值点对应的参数值ui之后，需要确定节点矢量U=(u0，u1，…，un)，最常采用的方法是平均值法。节点ui值为相邻的p-1个型值点对应的参数ui的平均值，有u0=…=up=0,um-p=…=um=1uj+p=1p∑i=jj+p-1u¯i，j=1,2,…,n-p （26）2.4.3 控制点求解在完成型值点的参数化和确定节点矢量之后，由式（27）易得到线性方程组（28），通过式（28）可以求得控制点Pi。qi=C(ui)=∑i=0nNi,p(ui)Pi （27）q0q1⋮qn-1qn=N0,p(u0)N1,p(u0)⋯Nn,p(u0)N0,p(u1)N1,p(u1)⋯Nn,p(u1)⋮⋮⋮N0,p(un)N1,p(un)⋯Nn,p(un)P0P1⋮Pn-1Pn （28）2.5　最小二乘法求解控制点一般情况下，可以通过型值点求出转子型线的控制点。但是，型值点增加后，控制点也会随之增加，这就会给后面的自动寻优带来一定的困难。因此，需要对控制点的个数进行优化。可采用最小二乘法减少控制点的个数[9-11]。最小二乘法是一种数学优化技术，通过最小化误差的平方和寻找数据的最佳函数。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。该方法在求解控制点时有着很大的优越性，具体步骤如下：1）型值点的参数化。2）确定节点矢量。3）求解k次基函数矩阵并将矩阵转置。4）计算曲线与型值点的误差。5）求解控制点。其中，对于控制点，可由下式求得(NTN)P=R （29）式中，N为K次基函数矩阵；P为控制点矩阵；R为误差矩阵。两边左乘(NTN)-1，可得P=(NTN)-1∙R （30）该式即为控制点矩阵的求解公式。3 型线分段及自动寻优3.1　型线分段如图3所示，本文以第三代典型转子型线中的复盛型线为例，其阴转子型线和阳转子型线都可以分为4段（图4、表1）。其中，阴转子型线由1段圆弧、2段圆弧包络线和1段椭圆弧包络线组成。而阳转子由1段圆弧包络线、2段圆弧和1段椭圆弧组成。先用一条4段5次的NURBS对其阴转子型线进行拟合，然后进行基于综合性能的优化设计。10.16578/j.issn.1004.2539.2023.11.008.F003图3复盛型线转子三维模型Fig. 3Three dimensional model of the rotor created by Fusheng profile10.16578/j.issn.1004.2539.2023.11.008.F004图4复盛阴转子型线分段Fig. 4Subsection of Fusheng profile on the female rotor10.16578/j.issn.1004.2539.2023.11.008.T001表1复盛型线分段Tab. 1Subsection of the Fusheng profile阴转子齿曲线阴转子曲线类型阳转子齿曲线阳转子曲线类型A2B2圆弧A1B1圆弧包络线B2C2圆弧包络线B1C1圆弧C2D2椭圆弧包络线C1D1椭圆弧D2E2圆弧包络线D1E1圆弧3.2　自动寻优要想进行自动寻优，就需要知道评判型线优劣的标准[12]。在双螺杆压缩机中，容积的密封性能是关键指标。其中，对容积的密封性能影响最大的是接触线长度L（mm）、泄漏三角形面积S（mm2）和面积利用系数M。下面分别介绍这3个性能参数。接触线[13]是指当双螺杆压缩机的阴、阳转子齿面啮合时，转子齿面发生接触，接触部分在齿面上形成的一条形状特殊的曲线。正是该接触线将基元容积的高压侧与低压侧分离。泄漏三角形面积[14]同样是一个非常重要的性能参数。通常情况下，双螺杆压缩机阴、阳转子接触线的最高点，如图5中的C点所示，无法到达机壳内壁面的交线（WW），因此，接触线最高点与阴、阳转子齿面以及机壳内壁面交线形成一个空间曲边三角形，如图5中的△ABC，称为泄漏三角形。其中，点A是阳转子与机壳内壁面交线的交点，点B是阴转子与机壳内壁面交线的交点，点C是阴阳转子接触线的最高点。10.16578/j.issn.1004.2539.2023.11.008.F005图5泄漏三角形空间位置示意图Fig. 5Spatial location diagram of the triangular leakage面积利用系数在性能评估中也占据着非常重要的地位。通常对转子型线进行优化设计时，减小阳转子齿间面积可能会导致阴转子齿间面积增大，因此，为更好地描述和量化优化设计前后阴阳转子齿间面积的变化，引入面积利用系数。面积利用系数计算式为M=z1(A01+A02)D12 （31）式中，z1为阳转子齿数；A01为阳转子齿间面积；A02为阴转子齿间面积；D1为阳转子齿顶圆直径。选择以上3个参数作为评判型线优劣的依据。其中，接触线长度和泄漏三角形面积与综合性能负相关，面积利用系数与综合性能正相关。目前已推出多级因素综合性能指标K=0.172 1L+0.725 9S/（0.102 0M），K值越小，型线的综合性能越好[15]。经验证，该综合性能指标能较好地反映容积的密封性能，精度略有提高，是可行可靠的。转子型线自动寻优的主要内容为：根据已知的型值点，结合最小二乘法对控制点进行求解；接着，基于综合性能指标K对控制点位置进行微调，以达到自动寻优的目的。其中，复盛型线原始控制点如图6~图9所示。10.16578/j.issn.1004.2539.2023.11.008.F006图6A2B2段数据点Fig. 6Data points on A2B210.16578/j.issn.1004.2539.2023.11.008.F007图7B2C2段数据点Fig. 7Data points on B2C210.16578/j.issn.1004.2539.2023.11.008.F008图8C2D2段数据点Fig. 8Data points on C2D210.16578/j.issn.1004.2539.2023.11.008.F009图9D2E2段数据点Fig. 9Data points on D2E2现在，分别向内外调整复盛型线的4段阴转子曲线，根据接触线长度、泄漏三角形面积以及面积利用系数来计算综合性能指标K，进而判断调整方式，具体如表2所示。10.16578/j.issn.1004.2539.2023.11.008.T002表2阴转子型线各段调整对综合性能指标的影响Tab. 2Influence on comprehensive preformance by adjusting the profile of the female rotor各段以及相应的调整方式K原始复盛型线1.003 2A2B2内1.002 7外1.008 6B2C2内1.007 7外0.962 9C2D2内0.998 0外1.003 9D2E2内1.003 2外1.003 1接着，向综合性能指标K值较小的方向进行调整，由表2可知，A2B2、B2C2、C2D2和D2E2的调整方向分别为向内、向外、向内和向外。从A2B2段选取21个数据点用以计算曲线的控制点，控制点的个数为7。其中，第一个控制点和最后一个控制点保持不变，对剩下的5个控制点进行调整，以达到优化的目的。在调整的过程中，需要用到控制变量法，即调整x坐标的时候，要保证y坐标不变；调整y坐标的时候，要保证x坐标不变；调整第二个控制点的时候，要保证第一个控制点保持不变。当各段曲线调整到最佳位置时，计算接触线长度、泄漏三角形面积、面积利用系数和综合性能指标，结果如表3所示。10.16578/j.issn.1004.2539.2023.11.008.T003表3阴转子型线各段调整到最优位置对综合性能指标的影响Tab. 3Influence on comprehensive performance by adjusting the subsection of the Fusheng profile of the female rotor to the optimal position调整线段接触线长度L/mm泄漏三角形面积S/mm2面积利用系数M综合性能指标K原始147.953 64.646 10.485 31.003 2A2B2147.990 04.508 40.487 30.908 0B2C2146.151 34.213 80.487 80.892 0C2D2146.482 34.205 70.495 90.970 5D2E2145.367 04.205 70.495 90.891 9其余各段型线的操作方法类似，不再赘述。由表3可知，相比于原始阴转子型线，优化后的型线在接触线长度、泄漏三角形面积、面积利用系数这3个性能参数上有所优化，具体对比如表4所示。10.16578/j.issn.1004.2539.2023.11.008.T004表4原始型线与优化型线性能参数对比Tab. 4Comparison of performance parameters of the original profile and the optimized profile接触线长度L/mm比例/%泄漏三角形面积S/mm2比例/%面积利用系数M比例/%综合性能指标K比例/%原始优化差值-1.75原始优化差值-9.48原始优化差值2.18原始优化差值-11.09147.953 6145.367 0-2.586 64.646 14.205 7-0.440 40.485 30.495 90.010 61.003 20.891 9-0.111 33.3　仿真分析现在，已经通过优化控制点位置得到了性能较优的阴转子型线，再通过相应的包络条件式，就可以生成阳转子型线。在UG软件中通过阵列和螺旋扫掠等功能，可以生成双螺杆的三维模型；然后，在ANSA软件中完成前处理；最后，导入到STAR-CCM+软件中进行相关设置并计算[16]。最重要的是设置双螺杆压缩机的进气压力和排气压力。其中，进气压力是一个标准大气压，排气压力为0.9 MPa。图10为型线寻优前后流场中的压力对比图。该图表明，在对原始复盛型线进行自动寻优后，双螺杆压缩机的压力由进气口向出气口递增时速率有所增加，超过了原始型线的双螺杆压缩机，且出口压力较大。这就表明型线自动寻优后的双螺杆压缩机的密封性能有所改善，从而验证了自动寻优理论的可靠性。图10压力分布对比Fig. 10Comparison of the pressure distribution10.16578/j.issn.1004.2539.2023.11.008.F10a1（a）型线寻优前的压力分布10.16578/j.issn.1004.2539.2023.11.008.F10a2（b）型线寻优后的压力分布 在对螺杆的实际加工过程中，还要综合考虑加工条件和运转性能，有时还要考虑噪声监测，从而对转子型线进行微调，进而加工出性能更加优越的螺杆转子。4 结论对双螺杆压缩机转子型线的正向设计展开研究，并基于NURBS曲线理论，对转子型线的控制点进行了求解。针对型值点法求解控制点会使控制点个数增加，从而给后续自动寻优带来困难的情况，提出了一种最小二乘法，可以有效减少控制点的个数，提高转子型线的设计效率。以复盛型线为例，基于综合性能指标对其进行自动寻优，并将原始型线和优化型线的性能参数进行对比，说明了具体的优化内容。最后，通过UG软件进行三维建模，利用STAR-CCM+软件对双螺杆进行CFD分析。结果表明，优化后的复盛型线与之前相比有更优的密封性能。