叶片曲面整体曲率光顺方法技术领域
本发明涉及微分几何以及叶轮叶片曲面造型方法，具体地，涉及一种叶片曲面整体曲率光顺方法。
背景技术
叶轮类零件在航空航天、能源动力等行业中均有广泛应用。叶片的三维造型是叶轮设计的难点，叶片曲面一般设计为自由曲面，型面结构非常复杂，而且叶片曲面需要满足几何连续性的要求，因为叶片曲面的光顺性对叶轮本身的效率有很大的影响。
目前叶片曲面的设计通常是先设计出叶片各流面上的流线，然后积叠流线得到叶片曲面。流线往往可以设计得很光顺，但由于没有对叶片叶高方向的光顺性进行控制，如曲率或曲率变化的控制，叶片沿叶高方向可能会出现皱褶，不光顺的情况。另外，叶片曲面的光顺性也会影响到后续的叶片数控加工工艺，如曲面不光顺会造成加工刀轨计算失败，或者生成的刀轨不光滑，造成加工过程中机床的振动。
发明内容
针对现有技术中的缺陷，本发明的目的是提供一种叶片曲面整体曲率光顺方法。
根据本发明提供的叶片曲面整体曲率光顺方法，包括如下步骤：
步骤S1：输入初始流线，插值初始流线得到叶片曲面；
步骤S2：计算叶片曲面在曲面采样点处沿流线方向和叶展方向的法曲率导数；
步骤S3：计算叶片曲面沿流线方向和叶展方向的法曲率导数平方总和E；
步骤S4：针对叶片曲面的点集{pk∈R3,1≤k≤m}中的每个点，计算每个点到原叶片曲面的距离d；其中，pk为点集中第k个点，R3为3维实向量空间，m为点集中点总个数；
步骤S5：若距离d大于给定的允许误差值dmax，则不调整该点，并在当前迭代次数未超过最大迭代次数Tmax时返回步骤S4以进行迭代计算；若距离d小于给定的允许误差值dmax，则利用黄金分割法搜索该点沿叶片圆周方向的最优调整角度θ，使叶片的法曲率导数平方总和E取最小值，并记为Enew；
若Enew＜E，则计算调整后该点的坐标，并替换原来的点，令E＝Enew，并在当前迭代次数未超过最大迭代次数Tmax时返回步骤S4以进行迭代计算；否则不调整该点，并在当前迭代次数未超过最大迭代次数Tmax时返回步骤S4以进行迭代计算。
优选地，所述步骤2包括如下步骤：
步骤S201：根据欧拉公式，叶片曲面在点p处沿切向x的法曲率计算公式κp(x)为
其中κ1表示曲面在点p处的最大主曲率，κ2表示曲面在点p处最小主曲率，表示切向x与主曲率κ1所在的主曲率方向的有向夹角；
步骤S202：由于叶片曲面是通过插值初始流线得到的，则认为叶片曲面的u参数线是叶片曲面的流线，v参数线方向是叶片的叶高方向；
步骤S203：设p1,p2两点是在叶片曲面的同一u参数线上，且两点的弧长距离为s，叶片曲面在点p1处沿流线方向的法曲率导数的计算公式为，
dκ p 1 ( t 1 ) d s = ( κ p 2 ( t 2 ) - κ p 1 ( t 1 ) ) s - - - ( 5 ) ]]>
其中t1表示叶片曲面在点p1处的流线方向，t2表示叶片曲面在点p2处的流线方向，表示叶片曲面在点p1处沿t1方向的法曲率，表示叶片曲面在点p2处沿t2方向的法曲率。
优选地，所述步骤3具体为：
根据沿叶片曲面流线方向和叶高方向给定的采样点数目，首先计算叶片在每点处沿流线方向和叶高方向的法曲率导数，进而计算叶片曲面沿流线方向和叶展方向的法曲率导数平方总和E，
E = Σ i = 1 n ( w 1 ( dκ p 1 ( t ) d s ) 2 + w 2 ( dκ p 1 ( l ) d s ) 2 ) - - - ( 6 ) ]]>
其中，表示叶片在点p1处沿流线方向的法曲率导数，表示叶片在点p1处沿叶高方向的法曲率导数；w1,w2分别表示叶片沿流线方向、叶高方向的法曲率导数平方总和权重；n为曲面采样点个数。
优选地，所述步骤S5中叶片曲面上点的调整方向被限定在圆周方向上，保证调整后的点与原来的点处于相同的流面上。
优选地，所述步骤5包括如下内容：
步骤S501：给定所述点的最大调整距离dmax，根据数据点的回转半径计算该点的最大调整角度θmax；在区间[-θmax,θmax]中，利用黄金分割法搜索最优调整角度θ，使得叶片曲面的法曲率导数平方总和E最小。
与现有技术相比，本发明具有如下的有益效果：
1、本发明解决了同时光顺优化叶片曲面沿流线方向和叶高方向的法曲率导数的问题，不仅能够使叶片曲面沿流线方向光顺，而且能够使叶片曲面沿叶高方向也光顺，使叶片曲面沿这两个方向的法曲率变化均匀，曲面变光滑；
2、本方法适用于轴流式和离心式叶轮叶片曲面的整体曲率光顺造型方法。
附图说明
通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：
图1为本发明的流程图；
图2为本发明中离心式叶轮叶片曲面在子午面上的流线模型示意图；
图3为本发明中计算叶片曲面在采样点处沿流线方向和叶高方向法曲率导数示意图；
图4为本发明中沿圆周方向调整叶片曲面数据点示意图；
图5为本发明中优化后的叶片曲面示意图。
图中：
1为在子午面上的叶片曲面流线；
2为叶片曲面上的采样点；
3为叶片曲面在同一流线上与采样点的弧长距离为s的点；
4为叶片曲面沿叶高方向与采样点的弧长距离为s的点。
具体实施方式
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进。这些都属于本发明的保护范围。
在本实施例中，如图1、图2、图3、图4以及图5所示，本发明利用整体优化法使叶片曲面沿流线方向和叶高方向的法曲率导数平方总和减少，使曲面变光顺。
输入初始流线，通过插值初始流线得到叶片曲面，根据欧拉公式计算叶片曲面在采样点p处沿流线方向和叶高方向的法曲率，
其中κ1,κ2表示叶片曲面在点p处的主曲率，表示切向x与主曲率κ1所在的主曲率方向的有向夹角。
计算叶片曲面在采样点处沿流线方向和叶高方向的法曲率导数。计算叶片曲面沿流线方向的法曲率导数公式如下，
dκ p 1 ( t 1 ) d s = ( κ p 2 ( t 2 ) - κ p 1 ( t 1 ) ) s - - - ( 2 ) ]]>
其中p1,p2两点是在叶片曲面的同一u参数线上，且两点的弧长距离为s，s是一个很小的数。t1表示叶片曲面在点p1处的流线方向，t2表示叶片曲面在点p2处的流线方向。计算叶片沿叶高方向的法曲率导数的方法与计算沿流线方向的类似。
根据给定的沿叶片曲面流线方向和叶高方向的采样点个数，计算叶片沿这两个方向的法曲率导数平方总和，
E = Σ i = 1 n ( w 1 ( dκ p 1 ( t ) d s ) 2 + w 2 ( dκ p 1 ( l ) d s ) 2 ) - - - ( 7 ) ]]>
其中，表示叶片在点p1处沿流线方向的法曲率导数，表示叶片在点p1处沿叶高方向的法曲率导数；w1,w2分别表示叶片沿流线方向和叶高方向的法曲率平方总和权重，可以根据两者的重要程度调整w1,w2的取值；
针对叶片曲面数据点集中的每一个点pi，如果其到原曲面的距离d。若d少于给定的允许最大距离值dmax，则利用黄金分割法搜索该点沿圆周方向的最优调整角度θ，使叶片的E取最小值；如果调整后的叶片曲面沿流线方向和叶高方向的法曲率导数平方总和比当前值小，则根据θ计算调整后该点的坐标，并替换原来的点，否则该点不作调整，对下一个点进行同样的操作。遍历完数据点集中的点后，进行下一次的迭代。如果当前迭代次数超过最大迭代次数，则程序结束。
本发明提供的叶片曲面整体曲率光顺方法，如图1所示。首先，输入初始流线，通过插值初始流线得到叶片曲面；计算叶片曲面沿流线方向和叶高方向的法曲率导数平方总和；对叶片曲面数据点集中的每个数据点，利用黄金分割法搜索该点沿圆周方向的最优调整角度，使叶片的法曲率导数平方总和取最小值；当遍历完数据点集中的点后，则进行下一次的迭代；如果迭代次数超过最大迭代次数，则程序结束，输出叶片曲面流线上的点。
在本实施例中，如图4所示的离心式叶轮叶片曲面数据点调整，类似的方法可以应用于轴流式叶轮叶片曲面，具体包括如下步骤：
步骤1：输入初始流线，通过插值初始流线得到叶片曲面；
步骤2：计算叶片曲面在曲面采样点处沿流线方向和叶展方向的法曲率导数；
κp(x)＝κ1cos2θ+κ2sin2θ(1)
其中κ1,κ2表示叶片曲面在点p处的主曲率，θ表示切向x与主曲率κ1所在的主曲率方向的有向夹角。
dκ p 1 ( t 1 ) d s = ( κ p 2 ( t 2 ) - κ p 1 ( t 1 ) ) s - - - ( 2 ) ]]>
其中p1,p2两点是在叶片曲面的同一u参数线上，且两点的弧长距离为s，s是一个很小的数。t1表示叶片曲面在点p1处的流线方向，t2表示叶片曲面在点p2处的流线方向。计算叶片曲面沿叶高方向的法曲率导数的方法与计算沿流线方向的类似。
步骤3：计算叶片曲面沿流线方向和叶展方向的法曲率导数平方总和E；
E = Σ i = 1 n ( w 1 ( dκ p 1 ( t ) d s ) 2 + w 2 ( dκ p 1 ( l ) d s ) 2 ) - - - ( 8 ) ]]>
其中，表示叶片在点p1处沿流线方向的曲率导数，表示叶片在点p1处沿叶高方向的曲率导数；w1,w2分别表示叶片曲面沿流线方向和叶高方向的法曲率导数平方总和权重，可以根据两者的重要程度调整w1,w2的取值；
步骤4：针对叶片曲面数据点集{pk∈R3,1≤k≤m}中的每个点，计算该点到原叶片曲面的距离d；
步骤5：若距离d大于给定的允许误差值dmax，则不调整该点，返回步骤4；若距离d小于给定的允许误差值dmax，则利用黄金分割法搜索该点沿叶片圆周方向的最优调整角度θ，使叶片的法曲率导数平方总和取最小值，并记为Enew；若Enew＜E，则计算调整后该点的坐标，并替换原来的点，同时令E＝Enew，返回步骤4；否则不调整该点，返回步骤4；
步骤6：如果当前迭代次数超过最大迭代次数Tmax，则结束，否则返回步骤4。
以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变形或修改，这并不影响本发明的实质内容。