本发明属于计算机辅助几何设计领域，涉及常微分方程、数值计算和曲线曲面造型等方面，是一种将新的表达形式用于造型设计之中的建模方法。
背景技术：
：制造业是国民经济的基础，在国家经济发展中占有至关重要的地位。复杂曲面工件高效高精度的数控加工已成为国家战略性装备与高新技术产业的迫切需求和制造技术的制高点。目前计算几何方法生成参数曲线曲面主要手段有：Bezier方法、NURBS方法等。Bezier方法要求具有插值和光滑拼接的苛刻条件，而NURBS方法，计算复杂，若选取权因子不适当，导致很坏的参数化，破坏曲面结构。具体而言，上述所提到的方法的处理对象都是数据散点，数据散点之上并没有考虑到曲面的动力学特性。例如，航空航天领域的压气机叶片和机翼等，部件处以复杂的空气流场中。这些流场对压气机叶片和机翼的设计起着至关重要的作用。理想情况下，这些部件表面的流场跟部件没有分离。可以运用部件表面流场来表示部件的形状，这样在极大程度上考虑到了曲面的动力学特性。然而，这样的问题带来了新的数据格式。每个散点处除了位置信息，还包含速度、压强、温度等。利用一阶常微分方程的做法可以较好地得到满足上述要求的叶片基元曲线表达形式，但由于常微分方程只能满足一个初始条件且刚性的特点，造型好的表达即使非常接近原曲面，其对误差和参数也极其敏感，而且叶片吸力面和压力面用一阶系统表达的基元曲线的末点无法重合，只能做到在一定误差的范围内接近，如果用于优化，一旦稍微调整系数矩阵的参数，曲线可能发生极大的变化，无法控制，从而失去优化的意义，只能作为造型之用。而二阶常微分方程的造型应用使得解决这些问题成为可能。技术实现要素：为了克服上述缺点，本发明从二阶常微分方程的具体形式出发，对压气机叶片形状进行高精度还原，提供了一种基于二阶常微分方程的压气机叶片吸力面基元曲线建模方法，并且获得了优化时较大的参数调节范围。本发明提供的基于二阶常微分方程的压气机叶片吸力面基元曲线建模方法，实现步骤为：步骤一：获取压气机叶片吸力面数据集，选用如下二阶线性非齐次常微分方程进行建模拟合。x··y··z··=A1xyz+A2x·y·z·+f(t)]]>其中：(x,y,z)是叶片在三维欧式空间中的坐标；A1和A2为系数矩阵；f(t)为关于时间t的多项式。设压气机叶片吸力面数据集表示为(X1,X2…,XN)，N为数据点总数，每个数据点为叶片在三维欧式空间中的一个坐标，其中第n个数据点Xn＝(xn,yn,zn)T，n＝1,2,…,N；采用如下差分形式进行数值拟合；xn+1+xn-1-2xnΔtn2yn+1+yn-1-2ynΔtn2zn+1+zn-1-2znΔtn2=A1A2xnynznxn+1-xnΔtnyn+1-ynΔtnzn+1-znΔtn+f1(tn)f2(tn)f3(tn)]]>其中，Δtn表示第n+1个数据点Xn+1和第n个数据点Xn间的步长，tn为第n个数据点处的累积步长，f1(tn)、f2(tn)和f3(tn)均为关于tn的多项式。步骤二：求取系数矩阵A和多项式f(t)。首先将压气机叶片吸力面数据集做如下排列：x3+x1-2x2Δt22...xn+1+xn-1-2xnΔtn2y3+y1-2y2Δt22...yn+1+yn-1-2ynΔtn2z3+z1-2z2Δt32...zn+1+zn-1-2znΔtn2=Ax2...xny2...ynz2...znx3-x2Δt3...xn+1-xnΔtny3-y2Δt2...yn+1-ynΔtnz3-z2Δt2...zn+1-znΔtn+f1(t2)...f1(tn)f2(t2)...f2(tn)f3(t2)...f3(tn)]]>其中，设矩阵矩阵矩阵则，进一步得到公式：D＝AM+F。考虑f(t)是多项式的情况，设系数矩阵为B，根据压气机叶片吸力面数据集，得到参数矩阵将F表示为F＝BT。则进一步地，获得如下形式：D=AM+F=AM+BT=ABMT]]>在上式的左右两边同时乘以其中是实对称矩阵，若是行满秩的，则是可逆的。等式两边同时乘以的逆，可以得到：AB=DMTT(MTMTT)-1]]>进而求得系数矩阵A。所求的为压气机叶片吸力面建立的模型可以表示为：其中X表示(x,y,z)T。步骤三：求基元曲线。根据边界条件的要求，将求取一个与所给数据首尾点重合的二阶曲线问题转化成为一个两点边值问题，并选取隐式欧拉法求解；设X(0)、X(1)表示压气机叶片吸力面数据集的起始点和终止点，求解如公式(6)所示问题得到基元曲线上的离散点。X··=AXX·+f(t)X(1)=xNyNzNX(0)=x1y1z1.]]>在求得基元曲线上的离散点后，进而可对基元曲线进行重构。本发明的优点与积极效果在于：本发明提出了一种以二阶常微分方程拟合压气机叶片吸力面数据散点的方法，解决了一阶系统无法同时保证首末端点插值条件的问题，获得了更大的优化操作空间，减小了对参数矩阵的敏感程度。本发明方法对于数据散点的首末点能够精确插值，这对于分片或分段拟合的数据散点的连接有着很大优势，可对压气机叶片形状进行高精度还原。附图说明图1是本发明的压气机叶片吸力面基元曲线建模方法的整体流程示意图；图2是叶片基元曲线(吸力面)拟合效果的正面示意图；图3是叶片基元曲线(吸力面)拟合效果的侧面示意图。具体实施方式下面将结合附图和实施例对本发明作进一步的详细说明。本发明一种基于二阶常微分方程的压气机叶片吸力面基元曲线建模方法，流程如图1所示，下面说明各个步骤。步骤一：获取压气机叶片吸力面数据集，对吸力面基元叶形数据进行建模拟合。(1)数据点的获得和预处理。本发明对压气机叶片初步设计得到的数据进行处理，将压气机叶片吸力面数据集导入。首先对数据进行预处理，包括：将得到的数据集统一尺度；通过适当的坐标变换，置于合适的坐标系下；去掉冗余或噪声数据等等。设预处理后得到的压气机叶片吸力面数据集表示为(X1,X2…,XN)，共有N个数据点，下面简称该数据集为原始数据集，原始数据集中的每个数据点为叶片在三维欧式空间中的一个坐标。原始数据集中的第n个数据点Xn＝(xn,yn,zn)T，n＝1,2,…,N。(2)对吸力面基元叶形数据进行建模拟合。本发明选用如下形式的二阶线性非齐次常微分方程进行拟合：x··y··z··=A1xyz+A2x·y·z·+f(t)---(1)]]>其中：(x,y,z)是叶片在三维欧式空间中的坐标；A1和A2为系数矩阵；f(t)为关于时间t的多项式，在本发明实施例中为三次多项式。在数值拟合时，可根据需要采用不同的差分格式。例如，原系统有如下差分形式：xn+1+xn-1-2xnΔtn2yn+1+yn-1-2ynΔtn2zn+1+zn-1-2znΔtn2=A1A2xnynznxn+1-xnΔtnyn+1-ynΔtnzn+1-znΔtn+f1(tn)f2(tn)f3(tn)---(2)]]>其中：{(xn,yn,zn)}为压气机叶片吸力面数据点集，n为数据点的序列号，n为正整数；设系数矩阵A＝(A1,A2)，只要求得A就可通过适当的初值条件计算出整条曲线。f1(tn)、f2(tn)和f3(tn)为第n个点处的非齐次项，均为关于tn的3次多项式。Δtn表示第n+1个数据点Xn+1和第n个数据点Xn间的步长，可以人为通过数据参数化方式给定，tn为第n个数据点处的累积步长，||Xn+1-Xn||表示数据点Xn+1和Xn之间的距离。步骤二，求取系数矩阵A和多项式f(t)。首先，将原始数据集做如下排列：x3+x1-2x2Δt22...xn+1+xn-1-2xnΔtn2y3+y1-2y2Δt22...yn+1+yn-1-2ynΔtn2z3+z1-2z2Δt22...zn+1+zn-1-2znΔtn2=Ax2...xny2...ynz2...znx3-x2Δt2...xn+1-xnΔtny3-y2Δt2...yn+1-ynΔtnz3-z2Δt2...zn+1-znΔtn+f1(t2)...f1(tn)f2(t2)...f2(tn)f3(t2)...f3(tn)---(3)]]>设矩阵矩阵矩阵则上面公式(3)可以写作下式：D＝AM+F(4)对于一般的f(t)，目前尚无求取的统一方法。本发明只考虑f(t)是多项式的情况。本发明实施例中以f(t)是一个三次多项式的例子说明。f(t)表示为如下：f(t)=b11b12b13t3+b21b22b23t2+b31b32b33t+b41b42b43---(5)]]>其中，(b11b12b13)T、(b21b22b23)T和(b31b32b33)T分别为系数向量，(b41b42b43)T为常数向量。将公式(5)写成矩阵的形式为：f(t)=b11b12b13b14b21b22b23b24b31b32b33b34t3t2t1---(6)]]>设矩阵设矩阵则进一步可得到：F＝BT。进一步地得到公式(4)有如下形式：D=AM+F=AM+BT=ABMT---(7)]]>在公式(7)的左右两边同时乘以其中是实对称矩阵，若是行满秩的，则是可逆的。等式两边同时乘以的逆，可以得到：AB=DMTT(MTMTT)-1---(8)]]>这样可求得A。所求系统为：X··=AXX·+f(t)=AXX·+Bt3t2t1---(9)]]>在具体求解时，不妨设A1＝0；执行下面过程：(2.1)计算步长Δtn和累积步长tn；(2.2)计算矩阵D和M；(2.3)计算矩阵T；(2.4)根据公式(8)计算(AB)；(2.5)重复上面过程(2.1)至(2.4)，获得吸力面15条基元曲线的3*7的矩阵为：(A1B1),(A2B2),...,(A15B15)。步骤三：求基元曲线。根据曲线上的边界条件，通过欧拉折线法还原出整条曲线。根据边界条件的要求，本发明希望得到一个与所给数据首尾点重合的二阶曲线，所以此问题转化成为一个两点边值问题，选取隐式欧拉法求解即可满足要求。即求解下列问题得到曲线上的离散点：X··=AXX·+f(t)X(1)=xNyNzNX(0)=x1y1z1---(10)]]>其中，X(0)、X(1)表示原始数据集的起始点和终止点。此处的X(0)、X(1)中的参数0和1表示的是累积的步长，对应的在起始点处累计步长为0和在终止点累计步长为1。具体来说，根据采用的具体形式不同，可以采用以下三类的公式：(1)A2＝0；此时，所求模型可表示为：利用下面公式求解：其中，I为单位矩阵，(2)A1＝0；此时所求模型表示为：利用下面公式求解：其中，(3)所求模型表示为：利用下面公式求解：上述三个公式(11)、(12)和(13)都是Sy＝b的形式。通过求解该线性方程组，得到中间离散点(X2,X3,…,XN-1)，也就完成了曲线重构。本发明实施例的一个具体求解过程如下：(3.1)以第一个点为二阶线性非齐次常微系统的初始点X(0)＝(x1,y1,z1)T，第71个点为末端点X(1)＝(x71,y71,z71)T。利用两点边值问题隐式解法，有其中，记为：SY＝b+f。计算得拟合后的基元曲线上的中间离散点(X2,X3,…,X70)。(3.2)重复上面(3.1)过程，得到吸力面15条基元曲线的拟合数据。本发明实施例中，将步骤三中获得的数据进行可视化，观察效果，如图2和图3所示。从图中可以看出以二阶常微分方法拟合的曲线精确插值数据散点的首末段，并且对于中间的数据散点的拟合也有很好的效果。当前第1页1 2 3