维修线性流量阀时的内筒孔设计.doc

1、全国第三届研究生数学建模竞赛题 目 维修线性流量阀时的内筒孔设计模型摘 要：油田固井机是石油重要开采设备，由于依赖进口整套设备价格昂贵，因此对固井机上受损的阀体部件进行修理具有重要的意义，可以带来很高的经济效益本文对维修线性流量阀是的内筒设计建立了三个模型，进行了详细的分析，并做了数值模拟，得出比较好的结果主要结果如下：（1）模型M-A(一般理论分析与简单模型）：建立了内筒孔为单参数曲线和双参数曲线的一般模型，从实际加工方面出发，对几种常见图形分别进行了讨论；（2）模型M-B(最小二乘意义下微元模型):通过微元法建立了一个最小二乘意义下的最优化模型，通过数值模拟，得到内筒孔是合适大小的矩形时，

2、线性化程度最高；（3）模型M-C（微分方程定量模型）：通过微元法建立了一个微分方程模型，内筒孔曲线与外筒孔圆的交点的横坐标与纵坐标所满足的微分方程组及其边界条件为：（1）初值条件： （2）通过数值计算与模拟，对问题一与问题二均给出了最佳内筒孔曲线设计方案，其大致形状为：参赛密码 （由组委会填写）关键词：线性流量阀；固井；数学模型；孔设计参赛队号 10414001 摘 要油田固井机是石油重要开采设备，由于依赖进口整套设备价格昂贵，因此对固井机上受损的阀体部件进行修理具有重要的意义，可以带来很高的经济效益本文对维修线性流量阀是的内筒设计建立了三个模型，进行了详细的分析，并做了数值模拟，得出比较好的

3、结果主要结果如下：（1）模型M-A(一般理论分析与简单模型）：建立了内筒孔为单参数曲线和双参数曲线的一般模型，从实际加工方面出发，对几种常见图形分别进行了讨论；（2）模型M-B(最小二乘意义下微元模型):通过微元法建立了一个最小二乘意义下的最优化模型，通过数值模拟，得到内筒孔是合适大小的矩形时，线性化程度最高；(3)模型M-C（微分方程定量模型）：通过微元法建立了一个微分方程模型，内筒孔曲线与外筒孔圆的交点的横坐标与纵坐标所满足的微分方程组及其边界条件为：（1）初值条件： （2）通过数值计算与模拟，对问题一与问题二均给出了最佳内筒孔曲线设计方案，其曲线大致形状为：关键词：线性流量阀；固井；数学

4、模型；孔设计1 问题的重述油田采油用的油井都是先用钻机钻几千米深的孔后，再利用固井机向四周的孔壁喷射水泥砂浆得到水泥井管后形成的.固井机上用来控制砂浆流量的阀是影响水泥井管质量的关键部件，但也会因磨损而损坏.目前我国还不能生产完整的阀体，固井机仍依赖进口.由于损坏的内筒已经被磨损得面目全非，根本无法测绘出原来的形状，因此维修时只能根据工作原理并结合阀的结构进行设计.根据仪表刻度可知控制流量的阀是一个线性阀，即阀体的旋转角度与砂浆流量成正比.在设计分析中假设砂浆的压力恒定，进而流量与“过流面积“成正比，因此阀体的旋转角度应该与“过流面积“成正比.一般来讲，控制流量的阀体为两个同心圆柱筒（两筒直径

5、大致相等）.外筒固定，它的侧面上有一个孔，形状为两个直径不等（相差至少3、4倍以上）的圆柱体的交线.内筒和外筒轴向之间没有相对运动，内筒可以自由转动.内筒的侧面上也有一个孔，它原来的形状未知（维修的任务就是设计内筒孔的形状），砂浆可以从两个孔的相交部分即“过流面积”流过.显然“过流面积”不能超过外筒孔的面积.现在数控机床比较普及，只要知道曲线的形状就可以在维修所需要的内筒上加工出合适的孔.当然从实际加工角度考虑，内筒孔的形状也不宜太复杂.第一个问题是讨论在上述阀体结构下，在“过流面积”从为零直到外筒孔面积的范围（简称“最大范围”）内，能否通过选择内筒孔形状实现“过流面积”与内筒旋转角度成严格的

6、线性关系.如果不能，要设计内筒孔的形状，在“最大范围”内，使“过流面积”与内筒旋转角近似成线性关系，同时在“最大范围”内，实际情况与严格线性关系的误差在某种意义下最小.第二个问题是考虑实际情况下，固井机向孔壁喷射水泥砂浆时经常采用的“过流面积”是在一个稍小的范围内，被称为主要工作区，它是 “最大范围”中的一段.因此，在维修固井机内筒时，比较令人满意的内筒孔形状应该使主要工作区中所对应的旋转角度的线性区间尽量长（至少达“最大范围”区间长度的75%以上），而且主要工作区的最大“过流面积”尽量大（至少要达到外筒孔面积的85%以上），并且使“过流面积”和内筒的旋转角度之间的“线性关系”尽量地好.请按此

7、要求设计内筒孔的形状. 2 问题分析 从已给的材料可以看出，内孔形状的设计有两个主要目标：（1）“过流面积”与内筒旋转角度的成线性关系程度越高越好;(2)尽可能得提高主要工作区间得最大“过流面积”.若将两个圆柱筒展开成平面，即为两个长方形，筒的转动可以转化为两个长方形的平动来考虑，那么目标（1）即等价于“过流面积”与圆柱筒平动距离的成线性关系程度越高越好.一般情况下，假设外筒固定不动，外筒孔曲线与内筒孔圆做相对平动，若其曲线交点始终不变，则移动的距离与“过流面积”是严格线性的.但是开始移动时刻其曲线交点不可能不变，因此一般不可能保持严格的线性，只能尽量的逼近线性,由于内外筒孔曲线相交初始阶段的

8、非线性化程度最高,因此设计内筒孔关键在于设计内筒孔曲线与外筒孔圆开始相交的那部分.3 模型假设 为使研究模型方便，本文做以下假设：1) 将圆柱筒展成平面,筒的转动转化为两个长方形平动,内筒旋转角与过流面积的之间的关系考虑为移动距离与过流面积的关系;2) 实际只要求过流面积从零到外筒孔面积范围内,因此假设最大过流面积并不一定要能达到外筒孔的面积;3) 内筒孔的形状简单易于加工,不能太复杂,是上下对称形状;4) 内筒孔的水平轴线与外筒孔的水平轴线完全重合,即两个孔的上下对称轴线要完全;5) 阀门保证过流面积从零到最大后不能继续再旋转.4 模型的建立、分析、求解与数值模拟4.1 模型M-A(一般理论

9、分析与简单模型)我们来考虑一般意义下的理论模型.外筒孔口的形状为一个圆，取（1,0）为坐标原点，圆的半径为1，则方程该圆的方程为： (1.1) 情形1 (适用于单参数曲线)设内筒孔口的前部形状（首先与外筒孔圆相交部分）为一满足方程的单参数封闭曲线，让该曲线向右平移得到一曲线族，该曲线族与外筒孔圆的交点的横坐标为，则两曲线所交面积为：，其中，是由方程所确定的函数且，显然与的关系并不一定是线性关系，现要通过选择适当的参数，使得尽量与成线性关系，即，从最小二乘法的思想出发，我们只要选择适当的参数使得极小即可，由多元函数极小的必要条件得：其中从而可以求解出得的值.情形2 (适用于双参数曲线)设内筒洞口

10、的形状为一满足方程的双参数封闭曲线，外筒洞口的形状为一满足方程的圆.让该曲线向右平移得到一曲线族，该曲线族与外筒洞口圆的交点的横坐标为，则两曲线所交面积为(即过流面积)：，其中，是由方程所确定的函数且，显然与的关系不一定是线性关系，现要通过选择适当的参数，使得尽量与成线性关系，即，由最小二乘的思想，我们只要选择适当的参数使得极小即可，由多元函数极小的必要条件得：其中，从而可以求解出得的值.实际上，如果不给出具体的曲线族，要对上述理论模型求解是非常困难的.因此我们通过合理的假设和附加一些约束条件来简化问题，比如规定好内筒孔与外筒孔首先相交的那部分的曲线形状（具体属于那类方程）,每次只要计算出每次

11、相交的点的横坐标，就可以得到交点处纵坐标，进而还可以求出过流面积与或者的关系.然后利用上面的理论模型，我们可以得到最优化条件下的或者.下面分别计算出内筒孔前部曲线分别为抛物线形、圆弧形、椭圆形、双曲线形的横坐标与过流面积.抛物线型，则，于是两曲线交点的横坐标为： (1.2)过流面积为:圆弧形，则，于是两曲线交点的横坐标为：, (1.3)过流面积为: ;椭圆弧型 则，于是两曲线交点的横坐标为： , (1.4)过流面积为:; 双曲线形，则，于是两曲线交点的横坐标为：,(1.5)过流面积为: 根据以上公式,我们可以进行数值模拟.有关此模型的几类形状的数值模拟结果见附录M-A.4.2 模型M-B（最小

12、二乘意义下微元模型）符号与变量说明: 外筒孔圆半径； : 微元单位步长； : 第个微元的面积； : 第个内外筒孔曲线交点纵坐标（上半部分）； : 第个内外筒孔曲线交点横坐标（上半部分）；模型分析将内筒壁和外筒壁展成平面后,转动角与过流面积之间的关系等价于相对移动距离与过流面积之间的关系.外筒壁上的孔是个圆,内筒前部边缘曲线的形状虽然未知,但不影响我们的讨论.如图(M-B-1)中所示,(其它曲线类型类似讨论,同时我们只考察圆上半部分,下半部分与上半部分对称).假设外筒孔圆的原心坐标为(0,0)，半径为，将圆形区域从左到右等成2n份，则每次前进步长为,其中用于分割的n充分大,交点集为从内筒孔以相对

13、外筒孔以为步长从左到右移动时,内筒孔前部边界与外筒孔圆边界的上半部分的交点.由于平动初始阶段过流面积与平动距离之间的关系受内筒孔右前部形状影响最大,因此我们只考察第二象限的1/4区间,来确定内筒孔右前部形状.该区间共有n+1个交点,起始交点为(-,0).由于每次前进的步长非常小,我们可以用梯形和三角形法计算每前进单位步长增加的近似面积为: (2.1)因此,显然我们要求面积与转动角度成线性关系,即要求该比值为常数,最可能的情况就是为常量.但内筒运动初始阶段不断在增加,因此这个比值不可能是常值,因此也不可能保持完全线性关系.理想情况下, 若要求过流面积从最小值到最大值能够保持严格线性关系,则每前进

14、,面积应该增加.模型建立通过以上分析,在最小二乘法意义下,线性度越高,即要求以下模型最优解:(MODEL M-B) (2.2)约束条件为: ; (2.3) (2.4). (2.5)当内筒相对外筒移动r时，各交点经过移动后最终停留位置的水平坐标为： (2.6)数值计算 这是一个非线性最优化问题,我们取,则,利用LINGO8.0求解,再利用MATALB6.1对数据进行处理（程序见附录M-B）,可大致判断出在最小二乘意义下内筒孔的形状.从结果分析可以看出内筒壁的孔的形状应该是一个矩形(忽略计算误差影响),对于, ,这个矩形的宽度(即高)为157.08,长度至少为200.显然,在这种意义下过流面积无法

15、达到外筒孔的最大面积.模型说明从一般常识也可以知道,要保持移动过程中做到近可能严格的线性化,内筒的孔中央的部分是一个矩形应该是合理的.因此虽然该模型最大过流面积无法达到外筒孔面积,但在最小二乘条件下,该模型也是与严格线性关系误差是最小的，从而给出了问题(1)的一个解答.但对于问题(2)，情况又有所不同.从模型M-B分析和数值模拟可以知道,内筒孔中间一段的形状为矩形时,过流面积能在一定范围内能够达到严格的线性关系.一般情况下,如图(M-B-3)所示, 内筒孔分为两部分,后部一段为矩形,前部为曲线与直线围成的部分(曲线尚未确定形状,也可能还是矩形).图(M-B-4)为内筒孔与外筒孔相交可能的效果图

16、(只考虑上半部分,下半部分对称).显然,从内筒孔与外筒孔相交于点到再次相交于点这个区间,过流面积与前进的步长是严格线性的. 设线段与的夹角为,现在考虑一种理想情况,按照模型（M-B）,内孔完全是个矩形的情形.按照要求,线性区间若要达到最大范围区间长度的75%,则的距离应该至少达到直径的75%,即要求: (2.7)此时图(M-B-2)中的长度为,的长度为,的长度约为,此时最大过流面积(即除去上下两个弓形部分圆内剩下的面积)与外筒孔面积比值应该满足: (2.8) 另一方面,若要主要工作区间的最大过流面积若要达到外筒孔的85%以上,则要求:,此时图(M-B-2)中的长度为,的长度为,的长度约为.此时

17、候的线性区间与最达范围的比值应该满足: .从以上分析可以看出,内筒孔完全设计为矩形无法同时满足以下两个条件:(1) 主要工作区间所对应线性区间至少达到最大范围的75%;(2) 主要工作区间的最大过流面积至少达到外筒孔的85%.因此,我们需要对内筒孔的前部曲线进行特殊设计建立更加精细模型M-C来满足问题(1)和问题(2)所有要求.4.3 模型M-C(微分方程定量模型)模型的分析与建立以(0,0)建立直角坐标系,则外筒孔圆的方程为.首先考虑这么一个问题:与圆交于点且斜率为k的直线,直线位置的微小变动对点横坐标的影响，即求点横坐标对直线在轴上的截距的微分（如图M-C-1）直线的方程为，向右平移得到与

18、圆方程联立得到与满足的隐函数关系由隐函数求导法则（3.1） (M-C-1) (M-C-2)假设内筒孔的上半部分是由函数决定的曲线当内筒孔曲线由左向右移动到与外筒孔圆相交于点时（外筒孔圆为），内筒孔顶点与外筒孔圆的左端点的距离为（如图M-C-2）过流面积为外筒孔圆与内筒孔曲线相交部分的面积此时内筒孔曲线的端点与交点的水平距离记为，的一个函数，因而交点的垂直位置为当内筒孔继续向右移动时，我们要考虑的是过流面积的局部变化特性为此，假设内筒孔曲线向右移动单位微元点水平移动到点，内筒孔曲线与外筒孔圆的交点变化为点内筒孔曲线与外筒孔圆的交点的变化记为内筒孔曲线在点的斜率为，外筒孔圆上点的坐标为由(3.1)

19、式可以知道由下式给出：从而（3.2）过流面积的局部变化微元为在外筒孔圆内部的内筒孔曲线扫过的部分与其在外筒孔圆上的交点变化所增加的曲边三角形面积之和（如图M-C-3）对于内筒孔曲线的微小平移，线段近似在过点的外筒孔圆的切线上由点坐标可以计算出点与点的垂直落差为， (M-C-3)而可以由三角形的面积代替，即从而（3.3）按照模型要求，过流面积应该与阀门转动角度成线性比例，即内筒孔曲线与外筒孔圆的相交部分的面积与交叉深度成线性比例增长因而要求须为常数或至少在某固定常数附近变化为此，我们令. （3.4）称出现在(3.4)中的常数为过流面积线性增长系数联立(3.2)式与(3.3)式得：简化上述方程组可

20、以得到：（3.5）考虑到在内筒孔曲线开始进入外筒孔圆时，内筒孔曲线与外筒孔圆相交于顶点处，即的微分方程组满足初值条件：（3.6）从而我们得到了在位置，内筒孔曲线与外筒孔圆的交点的横坐标与纵坐标所满足的微分方程组(3.5)及其边界条件（3.6）不同时刻点所描出的轨迹就能给出内筒孔曲线的形状及曲线数据数值求解与模拟我们利用Matlab(Version 6.5.1 Release 13 with Service Pack 1 )中对微分方程组(3.5)及其初值条件（3.6）运用Runge-Kutta(ode45)方法进行模拟（代码见附录M-C），结果数据分析如下（数值分析过程中，外筒孔圆默认为圆心在

21、原点上的单位圆） (M-C-4) (M-C-5)图（M-C-4）是为过流面积线性增长系数取1时，内筒孔曲线（左侧）及其刚开始进入外筒孔圆（右侧）的情形内筒孔曲线的形状呈对称流线型，顶部稍尖并且在内筒孔曲线进入外筒孔圆的过程中逐渐变钝呈带状我们需要指出的是这与我们前一节的分析是相容的在内筒孔刚开始要进入外筒孔时，由于接触到的外筒孔圆边缘呈竖直状态，内筒孔曲线需呈尖角状进入外筒孔，以满足交流面积增长当内筒孔与外筒孔圆交点处于外筒孔圆顶部（底部）时，圆的边缘处于水平状态，由交点改变引起的面积增长微元作用不明显，此时内筒孔曲线呈现秥滞型带状（函数数据见附表一）我们需要进一步分析的是上述内筒孔与外筒孔圆

22、逐步相交过程中，相交部分的面积与相交的深度（内筒孔曲线顶点与外筒孔圆左端点距离）的关系图(C-M-5)演示的是内筒孔曲线与外筒孔圆以步长为0.4，取值在0,2之间逐步重叠的过程每一步中实线围成的区域便是内筒孔曲线与外筒孔圆相交的部分(过流面积)回答问题一：我们将指出上述由微分方程组(3.5)及初值条件(3.6)给出的函数曲线基本满足维修要求的要求，过流面积与重叠深度呈线性比例关系我们首先还要指出的是，利用微分方程组(3.5)及初值条件(3.6)给出的数据，可以很方便地处理内筒孔曲线与外筒孔圆的位置关系附表中第一栏为ode45用自适应步长给出数值解时记录的内筒孔曲线的重叠深度，第二栏恰好是当内筒

23、孔曲线的重叠深度为时，内筒孔曲线与外筒孔圆的交点的横坐标，第三栏为纵坐标比如，当我们需要计算重叠深度时的过流面积，只需在第一栏中找到距最小的点和对应的第二栏数据在精度范围内便是重叠深度时内筒孔曲线与外筒孔圆的交点的横坐标当然，进一步对第二和第三栏用样条插值可以得到更高的精度下表(C-M-6)给出的是在取值于0,2之间，步长为0.1时，对附录数据中第二和第三栏用样条插值得到的数据进行数值积分后的过流面积关于改变的数据.t S(t)tS(t)0.10.20.30.40.50.60.70.80.91.00.0204472610.071044720.1492303980.2496210890.3681

24、141230.5013976940.646783275 0.8020711450.965465977 1.135494021.11.21.31.41.51.61.71.81.92.01.310957092 1.490865978 1.6744151471.860946167 2.049916409 2.240886354 2.433491554 2.627435758 2.822477744 3.01834542(M-C-6) 过流面积关于改变数据表结合图(M-C-5)，当即内筒孔顶点恰好到达外筒孔圆另一端时，阀到达最大位置此时内筒孔曲线完全进入外筒孔圆，过流面积达到最大最大过流面积达到外筒孔

25、圆面积的96.08图(M-C-7)是对表(M-C-6)中的数据进行最小二乘法拟合得到的拟合直线与表(M-C-6) 中的数据的对照关系图图(M-C-7)表明过流面积与重叠深度呈线性比例关系用数值积分可以得到表(M-C-6)数据的样条曲线与直线的偏差.我们称其为表(M-C-6) 中的数据的线性差当内筒孔体部继续进入外筒孔圆，过流面积继续增大一直到“最大范围”时，由图(M-C-5)内筒孔头部开始出即使这样，我们把上述的内筒孔尾部接上宽度为内筒孔最大宽度1.9634的矩形后，能很好地满足模型的要求回答问题二：观察图(M-C-7)我们不难发现，在内筒孔刚进入外筒孔圆时，过流面积呈现很严重的非线性性，而且

26、严重影响了全局的线性性质然而当在内筒孔的重叠深度不断增大时，过流面积呈规则的线性增长这是符合我们在开始所作的分析的在满足主要工作区大于75%的情况下，我们希望把主要工作区集中在阀门旋转的尾部 (M-C-7)图(M-C-8)反映的是将主要工作区放在重叠深度时，主要工作区过流面积和重叠深度的线性差随的变化规律显然只有当时才能达到最大范围区间长度的75此时过流面积最大值始终为最大过流面积达到外筒孔圆面积的96.08过流面积和重叠深度的线性差随主要工作区的减小而减小特别当主要工作区从2开始减小时，线性差迅速减小，但减小的趋势有减缓的趋向在实际中可以根据工作区与线性关系轻重关系取舍二者的矛盾因素 (M-

27、C-8)我们下面要考虑的问题从另一方也揭示了这种矛盾关系在上面求解微分方程组(3.5)及其初值条件(3.6)时，我们仅考虑的是将过流面积线性增长系数取1如果我们改变过流面积线性增长系数，会得到不同的结果下图(M-C-9)为过流面积线性增长系数取0.8时内筒孔曲线形状及其与外筒孔圆重叠的情况此时的最大过流面积，占外筒孔圆的76.9%但是此时对应的线性差降为0.12608 (M-C-9)图(M-C-10)表现的是最大过流面积和线性差对过流面积线性增长系数的变化趋势增长系数取值0.9,1；当增长系数超过1时，我们发现在重叠深度处，内筒孔会突出外筒孔圆，并且此时的微分方程组(3.5)往往变得不稳定最大

28、过流面积 (M-C-10)线性差我们进一步考察过流面积线性增长系数由图(M-C-7)可以发现过流面积在开始阶段增长较慢，因而我们可以在开始阶段增加过流面积线性增长系数，而在很小一段距离后恢复线性增长系数为正常值下面考虑在初始阶段重叠深度时，过流面积线性增长系数设置为1.25，之后恢复为1时，微分方程组(3.5)及初值条件(3.6)所给出的内筒孔曲线的数据及其分析图(C-M-11)演示的是内筒孔曲线与外筒孔圆以步长为0.4，取值在0,2之间逐步重叠的过程每一步中实线围成的区域对应内筒孔顶点所在位置的过流面积过流面积可以达到最大最大过流面积达到外筒孔圆面积的96.70图(M-C-12)是数据进行最

29、小二乘法拟合得到的拟合直线与过流面积数据的对照关系图,其中图(M-C-11) 图(M-C-12) 通过以上分析，我们改进的模型明显优于改经之前然而能不能以及怎样进一步通过对过流面积线性增长系数K的设置或把也K看成重叠深度t的函数，使得通过求解微分方程组(3.5)及其初值条件(3.6)得到更好结果但猜测，总有这样的对线性增长系数K的取法使得通过求解微分方程组(3.5)及其初值条件(3.6)得到的内筒孔能使得最大过流面积不断接近于外筒孔圆而不增加解的线性差5模型的评价与进一步讨论 在模型(M-A)中,我们建立了一般意义下的单参数和双参数的理论模型,并通过规定形状特征来求出关键函数,从而可以通过数值

30、模拟来确定最佳方案,最大的优点是:事先选取了常见的几何形状,特征固定,便于机床加工.缺点是,这样做无法保证方案一定是最优的.模型(M-B)通过微元法,建立微元模型,将问题转化为在最小二乘意义下的优化问题,通过数值模拟,得出了最小二乘意义下内筒孔最优方案是一个适当大小的矩形,主要缺点在于过流面积无法尽可能达到最大范围.模型(M-C)是一个建立在微分方程上更精细的模型,通过数值模拟我们可以发现,该模型很好的解决了主要工作区间所对应线性区间尽可能的长,主要工作区间的最大过流面积尽可能的大的要求,唯一不足的缺点是外筒孔曲线形状稍微有点复杂,但是随着数控机床的应用和不断发展,该不足可以得到弥补.当外筒孔

31、也磨损时，我们也可以考虑与外筒孔为圆的情形相似的办法在建立(3.1)式时，我们实际用到的是曲线上点的斜率当我们维修外筒孔发生了程度较轻的磨损的固井机时，可沿外筒孔边缘采集数据，进行插值拟合，得到相应点处的数值倒数数据，再按照建立微分方程组(3.5)式的过程建立微分方程组系统，用相似的数据分析可以求得对应内筒孔曲线的数据6参考文献1 姜启源，数学模型（第二版），高等教育出版社，1991. 2 叶其孝主编，大学生数学建模竞赛辅导教材（三），湖南教育出版社，1998. 3 朱道元，数学建模精品案例，东南大学出版社，1999. 4 雷功炎，数学模型讲义，北京大学出版社，1999. 5 萧树铁等，数学实

32、验，高等教育出版社，1999. 6 李尚志等，数学实验，高等教育出版社，1999.附 录M-A1.抛物线形 a0.20.40.60.81.01.2斜率k1.11541.42531.58031.65941.69931.7208截距b0.2510 0.3131 0.3343 0.3349 0.3261 0.31512.圆形a2832128256512斜率k1.75001.76661.76711.76711.76701.7670截距b0.2825 0.2224 0.2044 0.1998 0.1990 0.19863.椭圆(m,n)(2,3)(2,5)(2,9)(4,3)(4,5)(4,9)(7,3

33、)(7,5)(7,9)斜率k1.76461.76701.76711.75461.76611.76711.72551.76281.7670截距b-0.24000.2140 0.2031 0.2742 0.2288 0.2080 0.3112 0.2492 0.2151附 录M-B模型M-B相关代码:Lingo程序(求最小二乘意义下最优值)MODEL:SETS:dot/1.101/:y;ENDSETSDATA:ENDDATAmin=sum(dot(i)|i#NE#1:(y(i)+y(i-1)-3.1415926*100/2)2);for(dot(i):y(i)*y(i)0);y(1)=0;END%

34、(M-B)matlab程序(处理数据并绘出最终图形)data=0 0.00000014.10674 0.00000019.89975 0.00000024.31049 0.00000028.00000 0.00000031.22499 0.00000034.11744 0.00000036.75595 0.00000039.19184 0.00000041.46082 0.00000043.58899 0.00000045.59605 0.00000047.49737 0.00000049.30517 0.00000051.02940 0.00000052.67827 0.00000054.2

35、5864 0.00000055.77634 0.00000057.23635 0.00000058.64299 0.00000060.00000 0.00000061.31068 0.00000062.57795 0.00000063.80439 0.00000064.99231 0.00000066.14378 0.00000067.26069 0.00000068.34471 0.00000069.39741 0.00000070.42017 0.00000071.41428 0.00000072.38094 0.00000073.32121 0.00000074.23611 0.0000

36、0075.12656 0.00000075.99342 0.00000076.83749 0.00000077.65951 0.00000078.46018 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013

37、0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.

38、54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.000

39、00078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.00000078.54013 0.000000 78.54013 0.000000; y=data(:,1); y1=-y; x=-sqrt(10000-y.2); for i=1:101 x(i)=101-i+x(i); end theta=linspace(0,2*pi); plot(x,y,-,x,y1,-,(100)*cos(theta),100*sin(theta); axis equal;附 录M-C1 内筒孔曲线数据表25tl(t)F(t)=f(l(t)002.22E-162.51E-052.51E-055.02E-055.02E-055.02E-05