储油罐的变位识别与罐容表标定[优质内容]

1、2010高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A 我们的参赛报名号为（如果赛区设置

2、报名号的话）： A甲0701 所属学校（请填写完整的全名）： 青岛科技大学 参赛队员 (打印并签名) ：1. 唐坤 2. 蒋春林 3. 杨雪 指导教师或指导教师组负责人 (打印并签名)： 辛友明 日期： 2010 年 9 月 12 日赛区评阅编号（由赛区组委会评阅前进行编号）：2010高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：图表b储油罐的变位识别与罐容表标定摘要本文针对储油罐的变位识别与罐容表的标定问题，利

3、用投影积分法、近似替代法、多项式拟合以及误差修正函数建立了罐体变位后罐容表的标定模型。对于模型一，首先利用投影积分法分别求得无变位与纵向倾斜角度两状态下油位高度对应储油量的理论值；将此理论值与对应实际数据对比可得储油量误差，再分别对无变位时误差散点与两个状态误差的差值散点进行分析并拟合其曲线，由此便可确立的一次函数为修正函数并建立模型：最后通过符号积分进行模型求解。倾斜角度为时罐容表的标定值详见表1。对于模型二，首先利用投影积分法及倾斜球缺的油面近似替代法分别求得无变位与纵向倾斜角度横向偏转角度两状态下油位高度对应储油量的理论值；将无变位状态的理论值与其实验数据对比得储油量误差并拟合误差曲线，

4、由此建立含参数的修正函数并建立模型：然后利用计算机枚举搜索算法确定最小误差对应的值分别为，对应 的值分别为：，变位后罐容表的标定值详见表2。对于模型二的检验，可通过对比相同油位高度对应标定模型的理论值与对应实验数据，依据所得最大误差与总容量之比0.22%判断此模型较为准确。关键词： 投影积分 修正函数 拟合 计算机枚举搜索算法一、 问题重述通常加油站都有若干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表（即罐内油位高度与储油量的对应关系）进行实时计算，以得到罐内油位高度和储油量的变化情况。许多储油

5、罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化（以下称为变位），从而导致罐容表发生改变。按照有关规定，需要定期对罐容表进行重新标定。图1是一种典型的储油罐尺寸及形状示意图，其主体为圆柱体，两端为球冠体。图2是其罐体纵向倾斜变位的示意图，图3是罐体横向偏转变位的截面示意图。请用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题。 （1）为了掌握罐体变位后对罐容表的影响，利用如图4的小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为a=4.10的纵向变位两种情况做了实验，实验数据如附件1所示。请建立数学模型研究罐体变位后对罐容表的影响，并给出罐

6、体变位后油位高度间隔为1cm的罐容表标定值。（2）对于图1所示的实际储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度a和横向偏转角度b ）之间的一般关系。请利用罐体变位后在进/出油过程中的实际检测数据（附件2），根据你们所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。进一步利用附件2中的实际检测数据来分析检验你们模型的正确性与方法的可靠性。二、 问题分析由于测量时油位探针的下端固定于油罐底端，因此当罐体变位后所测罐内油位高度与储油量对应关系较变位前发生变化，即罐容表的标定值改变。为研究罐体变位后对罐容表的影响需分别

7、建立无变位和有变位两种状态下罐容表的理论标定模型，然后与实际测量值进行误差分析，对模型进行修正得到最终模型。2.1问题一的分析对于问题一，首先，考虑倾斜前后两种状态下的理论模型。对于无变位状态储油罐形状为两端平头的椭圆柱体，所以可直接将其简化为几何模型，将储油量转化为平卧椭圆柱体在某高度时的容积计算。此时油液所在立体为一规则的立体，可以通过简单的积分计算求得。对于变位倾斜角为状态，不考虑油浮子达到罐体顶端后对应储油量变化，并依据油位高度将倾斜后的几何模型分为四个部分计算。根据投影法思想，罐内油液上表面面积可用其在罐体侧面的投影面积计算，再通过微元法积分求解每一部分对应储油量，即得倾斜后罐容表的

8、理论值。其次，将所得油位高度对应储油量的理论值与对应实际数据对比得储油量误差，并对无变位时误差与两状态下误差之差分别进行分析并拟合其曲线，利用所得曲线分析倾斜角变化与相应的误差变化，并确立修正函数并进行模型修正，最后通过符号积分从而可得储油量和油位高度及变位参数的关系2.2问题二的分析对于问题二，要求建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度a和横向偏转角度b ）之间的一般关系。首先，考虑变位前后两种状态下的理论标定模型。变位前的理论模型可分为一个圆柱体与两个球缺体两部分并通过三重积分求解。对于横向偏转，通过垂直罐体底部的剖面图分析可知，不管是否有纵向倾

9、斜，横向偏转前测量油位高度与偏转后测量高度为一与横向偏转角b 有关的固定函数，因此变位后可先只考虑纵向倾斜再进行变量代换即可。而对于纵向倾斜时的测量高度与对应的理论储油量间的函数可类似模型一进行分段求解。由于罐体倾斜，液面与球缺的高度并不平行，但实际情况中倾角很小，因此可近似等效为球缺在垂直地面时液面为某高度的容积计算，可通过与的几何关系求的。球缺在垂直地面时液面为某高度的容积计算可通过简单的多重积分得出。罐身的容积与液高函数可同模型一利用投影法求的，由此便得变位前后的理论标定模型。其次，将无变位状态的理论值与其实验数据对比得储油量误差并拟合误差曲线，由此分析倾斜角变化与相应误差变化的关系并建

10、立含参数的修正函数对模型进行修正。然后利用计算机枚举搜索算法确定最小误差对应的以及的值，由此确定变为后罐容表的标定模型并利用实验所得数据隐形相关检验。三、 模型假设1）实际情况下储油罐变位角度不大。2）储油罐无严重变形，或有微小变形但可忽略不计。3）进出油量与油位高度的测量与时间点无关，或在短时间内可忽略其微小影响。4）所测油罐属性均为内测法测得，即可忽略油罐厚度对实验的影响。5）在试验的短期内温度、气压等没有剧烈变化。四、 符号说明符号说明储油量油位高度实际储油罐主体部分身长小椭圆型储油罐身长两次变位后的油位高度五、 模型的建立与求解5.1模型一测量时油位探针的下端固定于油罐底端，因此当罐体

11、变位后所测罐内油位高度与储油量对应关系较变位前发生变化即罐容表的标定值改变。为研究罐体变位后对罐容表的影响需分别建立无变位和变位倾斜角为两种状态下罐容表的标定模型，然后进行比较分析最终确立罐体变位后对罐容表的影响模型。5.1.1数据处理观察实验所得进出油数据，现不妨取无变位时两项数据作为研究对象。理论上某高度下的容积只与其油液高度有关，与其状态是进油还是出油无关，此处可利用定量分析进一步验证。 无变位出油所采集的数据为累加出油量而并没有给定出油前的总储油量，因此可利用数据中最大的油位高度对应于进口数据中累计进油量为，再加上罐内油量的初值便可得累加出油量为后油罐内剩余油量为，利用此方法可得累计储

12、油量对应剩余油量。由此可分别作出进出油过程油罐容量依据油位变化的曲线图一。图一：进出油过程油罐容量依据油位变化的曲线 由图一易观察出进出油过程两曲线重合，即某高度下的容积只与其油液高度有关，与其状态是进油还是出油无关，因此模型建立过程不考虑进出油的影响。5.1.2罐容表的理论标定模型 无变位状态 由图二可看出油罐形状为两端平头的椭圆柱体，由此将油罐储油量简化为几何图形的容积进行计算。图二：小椭圆型油罐形状及尺寸示意图1) 图二中阴影部分即为油罐油位高度为时对应储油的纵截面，其面积可利用定积分计算：依据图一建立的直角坐标可得油罐纵截面椭圆方程：由得的表达式，由此利用图一所示阴影部分的范围积分得其

13、面积为：2) 利用上述对其进行二重积分可得油罐油位高度为时对应储油量为：用定积分公式对求解得罐容表的理论标定模型：变位倾斜角为的状态 依据对应油位高度可将油罐分为图三所示五部分，其中罐体变位致使第5部分储油量范围内对应高度均为油位高度的最大值即此时油浮子到达油罐顶端，因此可将此罐体简化为只考虑前四个部分的几何模型。图三：油罐五部分及对应纵截面 第1部分1) 确定油位高度由图三易观察出储油量位于第1部分时油浮子位于罐体底端，即油位高度不随储油量增加而增加其值始终保持。2) 储油量的计算此时储油量与油位高度无关因此可直接计算第1部分的最大储油量。现利用微元法取油罐内厚度为的油层，其对应横截面为图四

14、中，在平面的投影为，则有：其中可利用椭圆方程进行定积分算出，则有：则有厚度为的油层体积为：而第一部分最大储油量油层厚度为，则有：图四：油层表面的投影 第2部分(1) 确定油位高度的变换范围 由图三易观察出第2部分油浮子开始随储油量增加而上升，可计算其最大储油量对应油位高度。由此确定油位高度的变化范围：(2) 储油量的计算此时储油量由两部分组成：一部分为第1部分的最大储油量；另一部分为罐内油的上表面到第2部分的下线间的油量，现对其计算：同样利用微元法取油罐第2部分内厚度为的油层，其对应横截面为图四，利用其在平面的投影为计算有： 首先利用定积分计算，利用图五(a)中确定积分上限为下线为，由此可得：

15、 则油罐内油层上表面到第2部分的下接线间的油量：并利用图五(a)中可确定上述积分上限为下限为，且有则可得此部分储油量：由此利用与确定储油量：图五：油层正面示意图 第3部分(1) 确定油位高度变化范围 第2部分所确定的变化上限为并且此处可作为油位高度变化的下限，而其下线可利用图三中的划分部分的相互关系计算得，由此确定油位高度的变化范围：(2) 储油量的计算此时储油量由两部分组成：一部分为第2部分的最大储油量；另一部分为罐内油的上表面到第3部分的下线间的油量，现对其计算：同样利用微元法取油罐第3部分内厚度为的油层，其对应横截面为图四，利用其在平面的投影计算有：首先利用定积分计算，如图五(b)中油面

16、对应映射点为，由可确定积分上限为下线为，由此可得：则油罐内油层上表面到第3部分的下线间的油量：，并利用图五(b)中,且可确定上述积分上限为下限为，则可得此部分储油量：由此利用与确定储油量： 第4部分1) 确定油位高度的变换化范围第3部分所确定的变化上限为且可作为此处油位高度变化的下限，而观察图三易知其上限为油位高度的最大值即，由此确立油位高度的变化范围：2) 确定储油量此时储油量由两部分组成：一部分为第3部分的最大储油量；另一部分为罐内油的上表面到第4部分的下线间的油量，现对其计算：同样利用微元法取油罐第4部分内厚度为的油层，其对应横截面为图四，利用其在平面的投影计算有：首先利用定积分计算，如

17、图五(c)中油面对应映射点为，由可确定积分上限为下限为，由此可得：则油罐内油层上表面到第4部分的下线间的油量：，并利用图五(c)中可确定上述积分上限为下限为，则可得此部分储油量：由此利用与确定储油量： 最终确立罐容表的理论标定模型5.1.3 修正模型的确立 理论值与实际值的对比 根据上述所得油罐储油量与油位高度的理论对应关系求解进油过程中实验所得油位高度的数据对应储油量的理论值，并分别在两种状态下将此理论数据与实验所得数据进行比较，并求得同一油位高度对应的理论储油量的值与实际储油量的值间的误差，如图六所示。图六：理论值与实际值的对比图图六中易看出每个状态的两条曲线并不完全重合，即油罐油位高度为

18、时对应储油量的理论值较实际值存在误差，且两个状态所产生的误差并不相同，但考虑到两种状态的测量环境和测量时间可能对两种误差引起的差别非常小，因此可猜测此误差的差距是由倾斜角度的不同直接引起的。因此下面着重分析倾斜角度对误差的影响。 误差曲线拟合 首先利用最小二乘法对无变位状态（=0）下所得误差散点进行拟合。拟合过程中通过对比几条拟合曲线的拟合优度最终确立无变位状态下的误差曲线为： 为研究倾斜角由0变到时上述误差产生的变化可将倾斜后所得误差与无倾斜时所得误差进行比较，并得到两种误差的差值。因此可对此差值的散点利用最小二乘法进行拟合，拟合过程中同样通过对比几条拟合曲线的拟合优度最终确立倾斜角改变前后

19、两误差的差值曲线： 此处即为由于倾斜角由0变到所导致的误差。图七：误差曲线拟合 确立修正模型 根据理论值与实际值对比过程中对影响误差的因素分析可假设此处误差只与倾斜角有关，则利用上述拟合函数建立储油量依据油位高度及倾斜角的误差模型：利用及对应的误差计算上式中系数，则可确立修正模型为：5.1.4罐体变位（倾斜角为）后对罐容表的影响模型利用上述所得倾斜角为时的罐容表的理论标定模型与修正模型可建立罐体变位后对罐容表的影响模型：5.1.5模型求解 利用积分符号对所得罐容表的标定模型计算油位高度间隔为的罐容表标定值见表1。表1：变为后油位高度间隔为的罐容表标定值油位高度（cm）罐容表标定值（L）油位高度

20、（cm）罐容表标定值(L)油位高度（cm）罐容表标定值(L)油位高度（cm）罐容表标定值(L)13.5331629.83611764.92913046.9226.5532660.76621807.50923087.49310.2433692.36631850.23933127.74414.3534724.59641893.08943167.65520.4635757.43651936.05953207.21627.5336790.87661979.13963246.39736.2337824.88672022.29973285.18846.1038859.44682065.52983323.5

21、5956.3939894.53692108.82993361.481070.1240930.13702152.161003398.941184.2441966.22712195.531013435.9212100.43421002.78722238.921023472.3813117.58431039.79732282.311033508.3014136.43441077.25742325.691043543.6615157.33451115.12752369.041053578.4116180.27461153.39762412.361063612.5417204.53471192.0577

22、2455.621073646.0018228.45481231.08782498.811083678.7719254.94491270.46792541.911093710.8020281.49501310.17802584.921103742.0521309.27511350.21812627.811113772.4722338.53521390.55822670.571123802.0123368.03531431.19832713.181133830.6224398.18541472.10842755.641143858.2225429.37551513.27852797.9111538

23、84.7326461.36561554.68862840.001163910.0627493.84571596.32872881.871173934.0528527.26581638.19882923.521183956.4529560.26591680.25892964.921193977.1230595.28601722.50903006.061203996.10 通过对比表1数据与实验数据发现在个别油位高度处仍对应储油量误差，但整体吻合较好，即所建模型准确性较高。5.2模型二5.2.1罐体无变位时罐容表的理论标定模型 球缺部分储油量的计算对于空间几何体（图八(a)中阴影部分）的体积可用三

24、重积分来计算，由图八(a)易看出几何体为球缺一部分，且由在坐标上的投影（图八(b)阴影部分）易知在轴上的变化范围为，再由在坐标上的投影（图八(c)阴影部分）只在轴上的变化范围为，在轴上变化范围为，由此可得：图八：球缺部分储油量示意图 确立罐容表理论标定模型圆柱体部分：可直接利用模型一中无变位状态油位高度为时储油量的计算方法，且此时有，则有：球缺部分：此时球缺部分油面高度即为油罐的油位高度，则有：得罐容表理论标定模型：5.2.2罐体变位后罐容表的理论标定模型 横向偏转（角度）现利用图九(a)分析横向偏转前后油位高度的关系。此正截面为以为半径为原点为的圆，且由于圆的特殊性使得偏转前后油面保持不变，

25、但偏转前后其油位探针所在直线由转变为，因此图九(a)中两条直线上的线段便可表示偏转前后油位高度。现依据高度分为三种情况讨论： 时此时有，其中，且，则可得： 时此时有，由可得： 时 此时有，其中，且，则可得： 且经验算只此时可作为三种情况的一般表达式。图九：横向偏转角度与纵向倾斜角度正截面图 纵向倾斜（角度）1）倾斜油罐圆柱体部分利用图九(b)将此油罐的圆柱体分为三个部分考虑，其对应储油量的求解方法与模型一中求解方法一致，因此可直接利用圆与椭圆的相似代换求得此圆柱体的储油量与油位高度的对应关系： 其中： 2）倾斜球缺部分 如图九(b)所示阴影部分分别为罐体倾斜后左、右球缺内储油量，可将其转化为无

26、变位时的储油量计算。由于倾斜角度非常小因此可近似认为图中分别为无变位时左、右两球缺达到的油位高度，因此可利用模型二中无变位时球缺储油量的计算公式进行计算。由图九(b)各线段分析知， ，代入上式得： 罐容表的理论标定模型的最终确立 倾斜角度时利用公式：总出油量=圆柱体储油量+球缺储油量可得 利用偏转角度前后油位高度关系式将上式中替换为可最终确立两次变位后储油量与油位高度的对应关系：5.2.3修正函数的确立 数据分析观察问题二中实际储油罐的采集数据，其中有出油量与显示油量容积两列，理论上分析任意采集时间对应两项数据之和应相同，但实际数据显示每个时间对应两数据之和均有差异。因此考虑油罐发生了倾斜或偏

27、转，而显示油量容积数据为油罐无变位时相应显示油高的对应值，由此可从采集数据中提取显示油高对应无变位时的油量容积以及对应发生倾斜或偏转时的储油量两类数据。 无变位状态下误差分析利用数据中显示油量容量作为无变位时显示油高对应的实际值，将其与上述所建罐容表的理论标定模型进行比较并计算相同油位高度对应储油量的理论值与实际值之差，将其作为误差。并利用最小二乘法对图十误差散点进行拟合。拟合过程中通过对比几条拟合曲线的拟合优度最终确立无变位状态下的误差曲线为：图十：储油量误差拟合曲线 确立修正函数忽略横向偏转角度对误差的影响，可依据模型一中修正模型的建立过程在此类似地建立修正函数，其中为上述建立的无变位状态

28、下的误差拟合曲线。由于此处偏转角度与倾斜角度均为未知量，则表达式可由模型一中变为后的误差散点图的大致走向确定为二次函数形式，即：。修正函数最终确立为：5.2.4建立模型并确定变位参数 利用上述所得罐容表理论标定模型与修正函数可确立罐内储油量与油位高度及变位参数间含待定系数的一般关系： 下面利用计算机枚举搜索算法求解误差最小时对应以及的取值，算法步骤：Step1 输入函数估算确定的搜索区间分别为：0，,0，-100，100, -100， 100, -100，100的取值为附录二给的实验测量值； Step2 当时,导入实验测量，得到相应值并分别与实验数据做差后平方再相加，得到并保存准确性指标；St

29、ep3 如果，转第二步并令，否则令并转第四步；Step4 如果，转第二步并令，否则并转第五步；Step5 如果，转第二步并令，否则并转第六步；Step6 如果，转第二步并令，否则并转第七步；Step7 如果，转第二步并令，否则转第八步；Step8 比较所有保存的，选择最小时所对应参数值作为最佳参数。 为确保所求数值的精确性可适当改变其搜索区间。通过编程实现上述算法并求解得：由此确立倾斜偏转后罐容表的标定模型：5.2.5模型求解及检验 模型求解利用积分符号对所得罐容表的标定模型计算油位高度间隔为的罐容表标定值见表2。表2：变为后油位高度间隔为的罐容表标定值油位高度（cm）罐容表标定值（L）油位高

30、度（cm）罐容表标定值(L)油位高度（cm）罐容表标定值(L)101302.3411021150.3621048751.38202053.2612023902.8722051248.09303279.1113026691.3423053636.22404862.3114029486.9624055871.47505691.1515032305.0725057912.05608771.6216035142.5126059718.787011003.1517037934.6627061254.818013304.2518040715.6928062407.839015867.5319043465.

31、8729063020.4110018483.0820046148.6230063581.75 模型检验对于模型二的检验可通过对比相同油位高度对应标定模型的理论值与对应实验数据之差得储油量误差，得图十一为其散点图。由图易观察出模型值与实验数据之差均保持在内，则绝对误差为最大误差与油罐的最大容量之比，计算得其值为0.22%，则说明模型二较准确。图十一：误差检验图六、 模型评价及推广6.1模型评价 模型优点： 1）在计算理论油液体积时运用了投影面积，等效近似等方法，使计算变得简单易行。 2）在求变位后罐容表的标定模型时，为抵消系统误差引入了与倾角有关的补偿函数对理论模型进行修正，是模型更加准确且具有

32、普适性。3）灵活运用了软件进行大规模的积分、数据拟合等运算。4）模型二建立过程中与模型一的修正标定函数联系分析，即简便了问题二得到求解，又使模型一的到了与推广。 模型缺点：1）模型将油罐在不同变位情况下的各种偏差（有由于倾角变化引起的也有由其他因素引起的），均转化为由倾角引起的情况，这就在某些情况下导致误差较大。2）在模型二中确定参数时，尚可大体估算区间，但对二次函数系数的区间估计则盲目性较大，只能尽可能扩大搜索区间，这就增大了计算机的搜索难度和搜索速度。6.2模型推广对模型进行适当改进后便可用于各种底面，各种倾角储油罐的容积计算与罐容表标定，并且可推广运用到其他各种涉及容器容积测定的行业。

33、参考文献1陈佰军 左振滨，编制卧式储油罐容量表的计算机模型，黑龙江八一农垦大学学报，第7卷 第2期：P55至P58页，1993。2战景林 王春平 王喜忠，倾斜椭平顶卧式罐容积的计算，中国计量，技术卷 检定、校准与测试：P73页，2006.4。3王郑耀，卧式加油灌剩余油料体积的计算，4国家技术监督局，中华人民共和国国家计量检定规程，JJG2261996。5欧阳光中 朱学炎，数学分析（第三版 下册），高等教育出版社，2008.4，附录程序一：%求变位前的实际、理论值x1=0.15902 0.17614 0.19259 0.2085 0.22393 0.23897 0.25366 0.26804 0

34、.28216 0.29603 0.30969 0.32315 0.33644 0.34957 0.36256 0.37542 0.38816 0.40079 0.41332 0.42576 0.43812 0.4504 0.46262 0.47478 0.48689 0.49895 0.51097 0.52295 0.5349 0.54682 0.55872 0.57061 0.58248 0.59435 0.60622 0.61809 0.62996 0.64185 0.65375 0.66567 0.67763 0.67854 0.69053 0.69082 0.70285 0.71491

35、 0.72703 0.73919 0.75142 0.7637 0.76416 0.77653 0.78899 0.80154 0.81419 0.82695 0.83983 0.85284 0.866 0.87932 0.89282 0.89284 0.90653 0.92045 0.93461 0.94905 0.9638 0.97891 0.99443 1.01043 1.02699 1.04425 1.06237 1.08159 1.10233 1.12532 1.15236 1.19349;y1=0.312 0.362 0.412 0.462 0.512 0.562 0.612 0.

36、662 0.712 0.762 0.812 0.862 0.912 0.962 1.012 1.062 1.112 1.162 1.212 1.262 1.312 1.362 1.412 1.462 1.512 1.562 1.612 1.662 1.712 1.762 1.812 1.862 1.912 1.962 2.012 2.062 2.112 2.162 2.212 2.262 2.312 2.31583 2.36583 2.36706 2.41706 2.46706 2.51706 2.56706 2.61706 2.66698 2.66883 2.71883 2.76883 2.

37、81883 2.86883 2.91883 2.96883 3.01883 3.06883 3.11883 3.16883 3.16891 3.21891 3.26891 3.31891 3.36891 3.41891 3.46891 3.51891 3.56891 3.61891 3.66891 3.71891 3.76891 3.81891 3.86891 3.91891 3.96891;a=0.6;b=0.89;d=2.45;y=(b*d/a)*(x1-a).*sqrt(x1*2*a-x1.2)+a*b*d*asin(x1/a-1)+(pi/2)*a*b*dt=y-y1plot(1000

38、*x1,1000*y1,r,1000*x1,1000*y,g,1000*x1,1000*t,*)legend(实际值,理论值,误差点值,Location,NorthWest);title(实际理论误差);程序二：%对倾斜角第一部分的积分syms x y a b t d1;a=0.6;b=0.89;d1=0.4;% t=(4.1/180)*pi; %t为纵向倾斜角 f1=2*b*sqrt(1-x2/a2);F1=int(f1,x,-a,y/cos(t)-a);V=int(F1/sin(t),y,0,d1*sin(t)%求的积分函数function p1=pp1(t)p1=-89/45000*(8

39、1*pi*cos(t)+16*(cos(t)2+3*cos(t)*sin(t)-1)/cos(t)2)(3/2)*cos(t)-108*asin(1/3/cos(t)*(2*sin(t)-3*cos(t)*sin(t)+162*asin(1/3/cos(t)*(2*sin(t)-3*cos(t)*cos(t)-108*(cos(t)2+3*cos(t)*sin(t)-1)/cos(t)2)(1/2)*cos(t)-54*pi*sin(t)/sin(t);%对倾斜角第二部分的积分syms x h a b d1 d2 t;a=0.6;b=0.89;d1=0.4;d2=2.05;% t=4.1/18

40、0*pi; %t为纵向倾斜角 f1=2*b*sqrt(1-x2/a2);F1=int(f1,x,-a,h+d1*tan(t)-a);V=int(F1*cos(t)/sin(t),h,0,h)%求的积分函数function p2=pp2(h,t)p2=-89/45000*(-162*asin(1/3/cos(t)*(2*sin(t)-3*cos(t)*cos(t)-16*(cos(t)2+3*cos(t)*sin(t)-1)/cos(t)2)(3/2)*cos(t)+108*asin(1/3/cos(t)*(2*sin(t)-3*cos(t)*sin(t)+108*(cos(t)2+3*cos(

41、t)*sin(t)-1)/cos(t)2)(1/2)*cos(t)-135*cos(t)*pi*h+2*cos(t)*(-25*h2-20*h*tan(t)+30*h-4*tan(t)2+12*tan(t)(3/2)-108*cos(t)*tan(t)*asin(5/3*h+2/3*tan(t)-1)-54*cos(t)*(-25*h2-20*h*tan(t)+30*h-4*tan(t)2+12*tan(t)(1/2)+162*cos(t)*asin(5/3*h+2/3*tan(t)-1)-270*cos(t)*asin(5/3*h+2/3*tan(t)-1)*h)/sin(t)+16.744

42、e-004; %对倾斜角第三部分的积分syms x h a b d1 d2 t;a=0.6;b=0.89;d1=0.4;d2=2.05;f1=2*b*sqrt(1-x2/a2);f2=2*b*sqrt(1-x2/a2);F1=int(f1,x,-a,h+d1*tan(t)-a);F2=int(f2,x,-a,h-d2*tan(t)-a);V=int(F1-F2)*cos(t)/sin(t),h,d2*tan(t),h)%求的积分函数function p3=pp3(h,t)p3=89/1440000*(2592*pi*cos(t)+343*(24*cos(t)*sin(t)-49+49*cos(

43、t)2)/cos(t)2)(3/2)*cos(t)-3024*(24*cos(t)*sin(t)-49+49*cos(t)2)/cos(t)2)(1/2)*cos(t)+21168*asin(1/12/cos(t)*(-49*sin(t)+12*cos(t)*sin(t)-5184*asin(1/12/cos(t)*(-49*sin(t)+12*cos(t)*cos(t)-5184*cos(t)*asin(5/3*h+2/3*tan(t)-1)-432*cos(t)*(-400*h2+1640*h*tan(t)+480*h-1681*tan(t)2-984*tan(t)(1/2)+5184*c

44、os(t)*asin(5/3*h-41/12*tan(t)-1)+1728*cos(t)*(-25*h2-20*h*tan(t)+30*h-4*tan(t)2+12*tan(t)(1/2)-64*cos(t)*(-25*h2-20*h*tan(t)+30*h-4*tan(t)2+12*tan(t)(3/2)+cos(t)*(-400*h2+1640*h*tan(t)+480*h-1681*tan(t)2-984*tan(t)(3/2)+8640*cos(t)*asin(5/3*h+2/3*tan(t)-1)*h+3456*cos(t)*tan(t)*asin(5/3*h+2/3*tan(t)-

45、1)+17712*cos(t)*tan(t)*asin(5/3*h-41/12*tan(t)-1)-8640*cos(t)*asin(5/3*h-41/12*tan(t)-1)*h)/sin(t)+0.1512553;%对倾斜角第四部分的积分syms x h a b d1 d2 t;a=0.6;b=0.89;d1=0.4;d2=2.05;% t=4.1/180*pi;f1=2*b*sqrt(1-x2/a2);F1=int(f1,x,-a,a-h+d2*tan(t);V=int(F1*cos(t)/sin(t),h,2*a-d1*tan(t),h)%求的积分函数function p4=pp4(h

46、,t)p4=-89/1440000*(5184*asin(1/12/cos(t)*(-12*cos(t)+49*sin(t)*cos(t)+5184*pi*cos(t)-1728*pi*sin(t)-3024*(24*cos(t)*sin(t)-49+49*cos(t)2)/cos(t)2)(1/2)*cos(t)+343*(24*cos(t)*sin(t)-49+49*cos(t)2)/cos(t)2)(3/2)*cos(t)-21168*asin(1/12/cos(t)*(-12*cos(t)+49*sin(t)*sin(t)-4320*cos(t)*pi*h+432*cos(t)*(48

47、0*h-984*tan(t)-400*h2+1640*h*tan(t)-1681*tan(t)2)(1/2)-5184*cos(t)*asin(-1+5/3*h-41/12*tan(t)-cos(t)*(480*h-984*tan(t)-400*h2+1640*h*tan(t)-1681*tan(t)2)(3/2)-17712*cos(t)*tan(t)*asin(-1+5/3*h-41/12*tan(t)+8640*cos(t)*asin(-1+5/3*h-41/12*tan(t)*h)/sin(t)+3.9589;程序三function y=th(h,t) %纵向倾斜角度%求分段函数的各个

48、理论值if h=2.05*tan(t) y=pp2(h,t);elseif h=2*0.6-0.4*tan(t) y=pp3(h,t);else y=pp4(h,t);end%模型一的修正函数 function y=xzh1(h,t)%变位之前对误差的拟合函数p3=-0.08403*h3+0.1506*h2+0.05822*h-0.001711; %变位前后对误差的差值的拟合函数q3=0.3925*h3-1.206*h2+0.9652*h-0.223; y=th(h,t)*1000-(p3+t/(4.1/180*pi)*q3)*1000;%作出修正函数的曲线clear all;clct=4.1

49、/180*pi;for i=1:120 a(i)=xzh1(0.01*i,t); b(i)=th(0.01*i,t)*1000;endplot(a,*);hold on;plot( b,*r);hold on; plot( a-b,*g);程序四function y=tou(h)%求取球缺的体积R=1.5;r=1.625;d=h/1000; %把油高一千等分t=0;y=0;for i=1:1000 u=sqrt(r2-(R-t)2); t=t+d; y=y+(r2-(R-t)2)*acos(r-1)/u)-(r-1)*sqrt(2*r-1-(R-t)2)*d;end%求取球身的表达式syms

50、x y z h R r LR=1.5;r=1.625;L=8;f1=1;f2=int(f1,x,-sqrt(R2-z2),sqrt(R2-z2);f3=int(f2,z,-R,h-R);f4=int(f3,y,-2,6)%求球体的总体积function V=zong(h)f4=9*pi+8*(-h*(h-3)(1/2)*h-12*(-h*(h-3)(1/2)+18*asin(2/3*h-1);V=f4+tou(h)*2%作出模型二的误差曲线a=xlsread(D:bb.xls,1,I2:J604);h=a(:,1);y1=a(:,2);R=1.5;r=1.625;d=h/1000;t=0;y=

51、0;u=sqrt(r2-(R-t)2);y=9*pi+8*(-h.*(h-3).(1/2).*h-12*(-h.*(h-3).(1/2)+18*asin(2/3*h-1)+(y+(r2-(R-t)2)*acos(r-1)/u)-(r-1)*sqrt(2*r-1-(R-t)2)*d)*2;plot(1000*h,1000*y1,r,1000*h,1000*y,g,1000*h,1000*(y1-y),*)legend(实际值 ,理论值,误差点值,Location,NorthWest);程序五% 作出区域一的总体积function V11=V11(h,t)syms h t r f H %t表示为纵向角度r=1.5;f=1/3;H=h+2*tan(t);V1=1/24*(-81*asin(1/3*(-4*sin(t)+3*cos(t)/cos(t)*cos(t)+16*2(1/2)*(3*cos(t)*sin(t)-2+2*cos(t)2)/cos(t)2)(3/2)*cos(t)+108*asin(1/3*(-4*sin(t)+3*cos(t)/cos(t)*sin(t)-54*2(1/2)*(3*cos(t