中国计量学院参赛论文-邓钧丞-蔡洪斌-于聪.doc

2010高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： A 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 中国计量学院 参赛队员 (打印并签名) ：1. 邓钧丞 2. 蔡洪斌 3. 于聪 指导教师或指导教师组负责人 (打印并签名)： 数模组 日期： 2010 年 9 月 13 日赛区评阅编号（由赛区组委会评阅前进行编号）：2010高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：21储油罐的变位识别与罐容表标定分析摘要本文主要对储油罐的变位识别与罐容表标定的问题进行了分析。对储油罐变位情况寻找出合理精确的标定罐容表的方法是很重要的。在问题（1）中，首先，我们根据试验中椭圆型储油罐的形状尺寸，用积分的方法建立了模型一，即未变位时的储油量与油位高度关系的模型。用模型一标定了一组罐容表作为未变位时的原罐容表，当发生倾斜角为的纵向变位后，我们用题目给出的实际检测数据，与原罐容表进行比较，发现存在很大的差异。然后我们对变位后的储油罐，同样用积分方法建立模型二，即倾斜角为的纵向变位后的储油量与油位高度关系的模型。考虑到实际储油方法中系统误差的存在，我们再对模型二进行改进，加入一个修正数，由实际数据拟合出修正数，得到了改进的模型二。再计算得到修正后的数据与实际数据之间的平均误差仅为0.11%，即验证了模型二精确度高。并用模型二计算标定出了变位后油位高度间隔为1cm的罐容表。在问题（2）中，我们首先根据储油罐倾斜后的情况，理论推出横向变位不改变容积的计算，只影响实际的油位高度。而纵向变位后，计算方法将发生改变，不再与未变位时相同。我们同第（1）问一样，用积分的方法建立模型三来计算储油量。实际模型中，储油罐的主体圆柱体内的容积计算方法类似于模型二的容积计算方法。对于两旁的球缺部分，我们根据液面重心近似法，将斜面近似转化成与重心在同一平面的水平面。此时的液面高度就是球缺部分液面的平均高度。之后，我们把顶板液面近似为一椭圆面并用积分法写出容积。然后，我们根据题目给出的出油量拟合出倾斜角度。拟合时，由于模型三的式子很复杂，无法直接拟合，我们用泰勒公式把式子展开化简为二次多项式，通过最小二乘拟合解得纵向倾斜角度为，横向倾斜角度为。将角度代入模型三，对罐体变位后油位高度间隔为10cm的罐容表进行标定，得到罐容表。接着，我们进一步利用附件2中的实际检测油高数据，用模型三求得储油量，算出各时的理论出油量，与实际的出油量比较，得到的理论值与实际值基本吻合，因此验证了我们的模型三是正确并且可靠的。最后，我们对模型的优缺点进行了评价，提出了模型改进的方向，并对模型进行了简单的应用与推广。关键词：变位识别；罐容表标定；泰勒展开；拟合一、问题的提出与重述1.1问题的提出通常加油站都有若干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表（即罐内油位高度与储油量的对应关系）进行实时计算，以得到罐内油位高度和储油量的变化情况。许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化（以下称为变位），从而导致罐容表发生改变。按照有关规定，需要定期对罐容表进行重新标定。附图1是一种典型的储油罐尺寸及形状示意图，其主体为圆柱体，两端为球冠体。附图2是其罐体纵向倾斜变位的示意图，附图3是罐体横向偏转变位的截面示意图。1.2 问题重述现用数学建模方法研究解决储油罐的变位识别与罐容表标定的如下两个问题：（1）为了掌握罐体变位后对罐容表的影响，利用如图4的小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为的纵向变位两种情况做了实验，实验数据如附表1所示。请建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。（2）对于图1所示的实际储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度a和横向偏转角度b ）之间的一般关系。请利用罐体变位后在进/出油过程中的实际检测数据（附表2），根据你们所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。进一步利用附件2中的实际检测数据来分析检验你们模型的正确性与方法的可靠性。二、问题分析本题要解决的是有关储油罐的变位识别与罐容表标定的问题。加油站的储油罐原先都有配套的“油位计量管理系统”，通常采用流量计和油位计来测量进/出油量与罐内油位高度等数据，来预先标定的罐容表。但许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生变位，从而使罐容表的值不再准确，需重新进行标定。经初步分析得：对于问题（1），要通过附图4的小椭圆型储油罐，来研究罐体变位后对罐容表的影响，并重新标定罐容表。可以先根据小椭圆型储油罐的形状和尺寸，运用高数知识推导出小椭圆型储油罐在未变位时的储油量与油位高度的关系式，以此标定出小椭圆型储油罐的罐容表，并以此表作为原先标准的罐容表。再对题目中给出的倾斜角为a=4.10的纵向变位的储油罐进行研究，由此时的位高所对应的储油量与罐容表对应的储油量进行对比做比较，看倾斜后的变化，再同样根据该储油罐的形状尺寸和倾斜度，运用高数知识推导出该储油罐此时的储油量与油位高度的关系式，根据关系式就可以给出罐体变位后油位高度间隔为1cm的罐容表标定值了。对于问题（2），要对附图1所示的实际储油罐，建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度a和横向偏转角度b ）之间的一般关系。当储油罐发生纵向倾斜角度a和横向偏转角度b的变位后，根据位置和形状尺寸，用数学积分的方法，在储油罐上取一微分容积，进行积分就可以求得油量关于探针测得的油位高度的函数关系式，即得实际储油罐的油容积模型。再利用附表2中罐体变位后在进/出油过程中的实际检测数据，根据所建立的数学模型确定变位后的纵向和横向的倾斜角度参数。这样就得到了此时位置下的储油量关于实测油位高的关系，然后根据该关系式给出罐体变位后油位高度间隔为10cm的罐容表标定值。最后再用附表2中实际的储油量数据来与此模型的计算值进行比较，根据对比结果可分析判断出此模型的正确性和方法可靠性。三、模型假设1、油位探针位置相对与油罐不发生变动。2、忽略油面张力导致的表面油接触罐壁时油面不平情况。3、忽略储油罐的壁厚，即储油罐的形状数据可视为油罐内部数据。四、符号及变量说明:为小椭圆型储油罐未变位的储油量(L)；：为小椭圆型储油罐变位后的储油量(L)；：为实际储油罐的储油量(L)；：为小椭圆型储油罐未变位的油位高(dm)；：为小椭圆型储油罐变位后的油位高(dm)；：为小椭圆型储油罐的罐长(dm)；：为实际储油罐的油位高度(dm)；：为实际储油罐发生横向倾斜后实际油面高(dm)；：为实际储油罐计算球缺顶部分容积的油面高(dm)；：为实际储油罐两侧球体顶点到圆柱体的距离(dm)；：为实际储油罐两侧球体的半径(dm)；：为实际储油罐主体圆柱体的直径(dm)。五、模型建立与求解5.1 问题（1）的模型建立与求解：加油站的储油罐都有对应的罐容表，此罐容表是在储油罐水平无变位时测出的，但使用一段时间后，由于地基变形等原因，使罐体的位置会发生变位。为了研究罐体变位后对罐容表的影响，在这一问，我们利用小椭圆型储油罐，建立了分别在未变位时和变位后的两个罐容表标定的模型进行研究。5.1.1模型一（未变位）：A模型一的建立小椭圆型储油罐为两端平头的椭圆柱体，储油罐中油的体积是油面高度的函数，可表示如下： （1）式中：为储油体积（L），为油位高度（mm）。a-ab-b油面截面dy（x,y）h图5 椭圆罐侧剖图由图5，根据积分概念，体积元素为： （2）式中：为截面面积，是关于的函数，是截面与椭圆相交时椭圆的纵坐标。的表达式如下： (3)式中：为截面与椭圆相交时椭圆的横坐标，为储油罐的纵向从长度。 椭圆方程为： （4）式中：为短半长轴，为长半长轴。再有与油面高度的关系为: (5)由以上（1）.（5）5个式子可以解得储油罐中油的体积和油面高度的函数模型如下：模型一 B模型一的求解与结果检验分析运用MATLAB软件，分别带入附表1中未变位的进油与出油时测得的油位高度值求得对应的储油量。再分别由附表1中未变位的进油与出油量算得相应各油位高度时的实际储油量，做出对应油位高度的理论算得的储油量与实际的储油量的分布图如下图6所示：图6 理论算得的储油量与实际的储油量的分布图由结果可以看出：（1）进油与出油在两种情况中都几乎重合，可以认为进油与出油时，储油量与油位高度的关系相同；（2）相同油位高度时，理论计算出的储油量比实际的储油量稍大，且随高度的增加，差值增大。C模型一改进与改进后检验根据上一步的结果分析，理论值与实际值之间存在一定的误差，且随油位的增加而增大，可以认为这个误差包含公式精度、随油位的增加压强增大导致的体积减小、储油罐的不规则、探针和进出油管的体积等因素。对模型一进行改进，因为差值随油位的增加大致成线性增加，我们给模型一加一个修正数，如下： （6）式中：为修正数，是关于油位的函数，如下： （7）式中：c、d为参数。用实际值和理论计算出的值进行拟合，用MATLAB解得： (8)改进后的模型一：加上修正数后，进油与出油时的理论值与实际值分布如图7所示：图7 油与出油时的理论值与实际值分布图 由图7可以得到，改进后的模型求解的值与实际值几乎吻合，且平均误差为0.57%，认为该模型较为精确。用此模型一做出小椭圆型储油罐未变位时的罐容表如附表3所示。5.1.2模型二（变位后）：当小椭圆型储油罐发生倾斜角为的纵向变位后，观察附表1中倾斜变位后进出油的各油位高时的储油量，与附表3中未变位的原罐容表的对应相同油位时的储油量，可以发现，倾斜变位后的油位高所对应的实际储油量与原罐容表所对应的储油量发生了很大变化，原罐容表不再适用于倾斜变位后的储油罐了，需重新标定倾斜变位后储油罐的罐容表。A模型二的建立小椭圆型储油罐发生倾斜角为的纵向变位，我们取它的正视图建立坐标系，以储油罐的正视图的边位坐标轴，并由倾斜的角度，根据油面的不同情况将油面在储油罐中分为5种情形进行分析求解，如图8所示：A(0,0)B(4,0)C(24.5，0)D（24.5，11.2438）E（24.5，11.5305）F（24.5，12）G(4,12)H(0,12)I（0，1.7562）G（0，0.2867）K(4,1.469)L(4,11.7133)油面图8 油罐正面剖视图假设油位探针所测油位高为，单位为m,储油量为，单位为L。对5种情况分别计算：第1种情况：当油面刚达到或未达到BG时，此时刚达到或未达到油位探针的最低端B点，油位探针所显示的油位高为0m，此时有最大时为油面与BG重合时的储油量，由于油面在BG以下时，探针不能测出，这里只要知道最大的油量即可。此时有： （9）第2种情况：当油面在BG与CI之间时，油面在探针BK段，所测油位高的取值范围为,如下图9所示：A(0,0)B(4,0)C(20.5，0)DEFG(4,12)H(0,12)IGK(4,1.469)L(4,11.713)油面dy图9 第2中情况时的油罐正面剖视图根据图有： （10）第3种情况：当油面在CI与DH之间时，油面在探针KL段，所测油位高的取值范围为,如下图10所示：A(0,0)B(4,0)C(20.5，0)DEFG(4,12)H(0,12)IGK(4,1.469)L油面MN图10第3中情况时的油罐正面剖视图由图可以以MN为界将容积分为两部分，即： (11)其中：为上面部分容积，为下面部分容积。计算与第2种情况一样，的计算同模型一的计算方法一样。则有： （12） （13）第4种情况：当油面在DH与EG之间时，油面在探针LG段，所测油位高的取值范围为，如下图11所示：A(0,0)B(4,0)C(20.5，0)DEFG(4,12)H(0,12)IGK(4,1.469)L(4,11.713)油面PQ图11第4中情况时的油罐正面剖视图由图可以以PQ为界将容积分为两部分，即： (14)其中：为上面部分容积，为下面部分容积。计算与第2种情况一样，的计算同模型一的计算方法一样。则有： （15） （16）第5种情况：当油面刚达到或超过EG后，此时油面刚达到或已超过油位探针的测量范围，油位探针所显示的油位高为1.2m，此时储油量为： （17）通过以上分析可得储油量的关系如下模型二：第1种情况： 第2种情况：第3种情况： 第4种情况： 第5种情况： B模型二的求解与结果检验分析 由模型一的结果和分析知道，进油与出油的油位高度与储油量间的关系没有变化，这里我们就只考虑用进油的数据进行求解。运用MATLAB软件，带入附表1中倾斜变位后的进油时测得的油位高度值求得对应的储油量。再由附表1中倾斜变位后的进油量算得相应各油位高度时的实际储油量，做出对应油位高度的理论算得的储油量与实际的储油量的分布图如下图12所示：图12 理论算得的储油量与实际的储油量的分布图 由结果可以看出：理论计算出的值与实际值之间存在较大的误差。C模型二改进与改进后检验根据上一步的结果分析，理论值与实际值之间存在一定的误差，即理论进行计算时，相比实际情况还存在一些误差。对模型一进行改进，加一个修正数，即： （18）式中：为修正数。计算出实际值与理论计算出的值之间的误差，做出误差随油位高度的分布，并拟合出得到的表达式如下： （19）则改进后的模型二为： （20） 运用改进后的模型二进行求解得到一组理论计算值，再与实际值做分布进行对比，如图13所示：图13 理论计算值与实际值分布图 由图可以得到，改进后的模型求解的值与实际值几乎吻合，且平均误差为0.11%，认为该模型精确。用此模型二对小椭圆型储油罐倾斜a=4.10变位后的罐容表进行重新标定，重新标定后罐容表如表4所示：表4油位高（dm）储油量（L）油位高（dm）储油量（L）油位高（dm）储油量（L）00,1.67154965.6582661.40.13.5254.110058.12703.50.26.25974.21044.68.22745.50.39.96874.31084.58.32787.20.414.7534.41124.88.42828.70.520.6864.51165.38.528700.627.8484.61206.28.62911.10.736.3134.71247.28.72951.80.846.1394.81288.68.82992.30.957.3894.91330.18.93032.5170.12251371.993072.41.184.395.11413.89.131121.2100.255.214569.23151.21.3117.745.31498.39.33190.11.4136.925.41540.89.43228.61.5157.825.51583.59.53266.71.6180.255.61626.39.63304.41.7203.995.71669.29.73341.71.8228.95.81712.29.83378.51.9254.885.91755.39.93414.92281.8561798.5103450.72.1309.756.11841.810.13486.12.2338.546.21885.110.23520.92.3368.146.31928.510.33555.12.4398.526.41971.910.43588.82.5429.656.52015.410.53621.82.6461.486.62058.810.63654.22.7493.996.72102.310.73685.92.8527.146.82145.710.83716.92.9560.96.92189.110.93747.23595.2472232.5113776.63.1630.147.12275.811.13805.33.2665.587.22319.111.238333.3701.527.32362.311.33859.83.4737.957.42405.411.43885.63.5774.857.52448.411.53910.33.6812.27.62491.311.63933.93.7849.977.7253411.739563.8888.157.82576.611.83976.63.9926.717.92619.111.93995.5124012.7,4110.15.2 问题（2）的模型建立与求解：根据问题（1）的分析，当储油罐发生变位倾斜后，探针所测的油位高度所对应的实际储油量与对应原未变位时的罐容表所示的储油量之间有很大的差别，原罐容表不再适用倾斜变位后的，需重新进行罐容表的标定。对于附图1所示的实际储油罐，我们建立罐体的位置为一般情况下，即纵向倾斜角度a和横向偏转角度b变位，测量的油位高度与储油量的函数关系模型三，再由此模型来标定倾斜后的罐容表。A模型三的建立附图1所示的实际储油罐，其主体为圆柱体，两端为球冠体。由附图2可以看出，纵向倾斜后，不仅对测量油面高度产生影响，还对容积的计算产生影响。由附图3可以看出，横向倾斜时，由于实际的储油罐，其主体为圆柱体，两端为球冠体，所以不对容积的计算产生影响，只对测量油面高度产生影响。横向偏转角度b变位，纵向倾斜角度a变位，侧面截面图如下图14所示：oHH1Rb油面探针图14 侧面截面图则有： （21）式中：实际油面高，为测量的油面高，为储油罐圆柱体的半径。正面截面图如下图15所示：h1油面探针VaVbVc图15 正面截面图则有模型三，即总油量为： （22）式中：为左边球内的油量，为中间在圆柱体内的油量，为右边球内的油量。对于的计算，与模型二的类似，只是这里是圆柱体，模型二是椭圆体，也可分为5种情况，这里我们不再细分列出，而用一个总式子表示，即： （23）对与和的计算，由于罐体倾斜，使顶型部分的液面也相对为一个倾斜面，若采用水平状态计算其部分容积，此时的液面高度就应该是液面的平均高度。我们把顶板液面近似为一半椭圆面，根据文献1，得其重心在罐体轴先方向的距离为（为短半长轴）。如图16所示；若高端液面高为时，顶型体内液面半椭圆的重心为G，则： (24)就可用OE液面来近似代替OF液面来计算液面以下顶型部分的油量。油量高为： （25）H1CGFEOIAH2图16 正面截面图球顶部分容积的计算，如图17所示：D/2edexlROH2dee图17 球缺顶顶型部分容积计算原理图取球缺体上一微分容积，即： （26）因为 (27)所以 （28）球缺体在高度下的部分容积： （29）运用泰勒公式展开得：其中：,，。 根据式子（21）、（13）（29），带入（22）式，即得模型三。B模型三的参数求解与罐容表的标定由于模型三的式子很复杂且很长，用计算机直接进行拟合无法算出结果。我们先对模型三的式子用泰勒公式展开化简，得到一个相对原来较简化的式子：其中：然后根据附表2中的出油量和油位高度两组数据，对模型三中的两个倾斜角度进行拟合，即：通过MATLAB拟合得到：，。将角度代入模型三中，再对罐体变位后油位高度间隔为10cm的罐容表进行标定，标定的罐容表如表5所示：表5油位高度（dm）储油量（L）油位高度（dm）储油量（L）油位高度（dm）储油量（L）00101438820414611143111689921441372438.71219495224674831065.91322159234927942275.41424876245171453804.51527631255403465577.51630409265622077551.11733196275825189693.618359762860101911980193873629617433063241C模型三的正确性和可靠性的分析根据附表2中显示油高的数据，代入模型三，计算出对应的储油量，再由储油量算出油高差值对应的理论出油量，做出与附表2中出油量的分布图，如图18所示：图18 实际出油量与理论出油量分布图 由图18可以看出，理论算出的出油量与实际的出油量基本吻合，既可以验证模型三是正确的，并且可靠。六、模型检验 此部分内容见第五步模型建立与求解中：5.1.1B、5.1.1C、5.1.2B、5.1.2C和5.2C步骤。七、模型评价与改进模型的优点： 1、使用积分的方法求得的容积数据精确；2、计算方法简单，易于理解；3、通过泰勒展开使得一个看似无比巨大的积分化为具有相当精度的二次展开式。模型的不足之处：对误差的处理不够。可行的改进方向：在实际生活生产中，我们不能随时高精度测量计算误差，有时我们不得不进行误差较大但可行度较高的近似算法，此处我们推荐泰勒的一次线性展开式，通过计算可得： 带入实际数据，发现在已知，情况下，实测液高可以很好反映储蓄罐中的液容量。八、模型的应用与推广本文提出的模型，是使用积分的方法来求体积，结果精确，对类似的储油罐的容量标定都可以应用。特别是一般情况的模型三，可以应用到很多场合。 对于一般其他形状的储油罐，或倾角方向改变，也可以用本文提出的方法类似进行计算。参考文献：1廉育英，容量计量技术M，北京：中国计量出版社，2006。2刘慧颖，MATLAB R2006a基础教程M，北京：清华大学出版社，2007。3同济大学数学系，高等数学M，北京：高等教育出版社，2007。4曾强鑫，油品计量员培训教程M,北京：中国石化出版社，2005。附录：程序：模型一部分程序：h=load(h.txt);h=h./100;v=;for i=1:78 (h)2*8.9/6*24.5*(62-(h-6)2)(0.5); fun=(h)2*8.9/6*24.5*(62-(h-6)2)(0.5); v(i)=quadv(fun,0,h(i);endvh1=load(h1.txt);h1=h1./100;v3=;for i=1:74 (h1)2*8.9/6*24.5*(62-(h1(i)-6)2)(0.5); fun=(h1(i)2*8.9/6*24.5*(62-(h1(i)-6)2)(0.5); v3(i)=quadv(fun,0,h1(i);endv3v2=load(v2.txt);v1=load(v1.txt);v1=v1+262;plot(h,v1,h,v,，h1,v2,*,h1,v3,*)V=v-v1; H=h; V=V; H=ones(78,1) H; b,bint=regress(V,H) V1=v+b(1).*h+b(2); plot(h,v1,h,V1); V2=v3+b(1).*h1+b(2); plot(h1,v2,h1,V2)plot(h1,v2,h1,V2,h,v1,h,V1)模型二部分程序：for i=1:15y2(i)=0.01*(i-1);fun=(y)2*0.89*sqrt(1-(y-0.6)2/0.36)*(-y/tan(2*pi*4.1/360)+0.4+y2(i)/tan(2*pi*4.1/360);v0(i)=quadv(fun,0,y2(i)+0.4*tan(2*pi*4.1/360);endv0for i=1:103 y2(i)=0.15+0.01*(i-1); fun1=(x)4.9*0.89*sqrt(1-(x-0.6).2/0.36); fun2=(y)2*0.89*sqrt(1-(y-0.6)2/0.36)*(2.45-(y-y2(i)+2.05*tan(2*pi*4.1/360)/tan(2*pi*4.1/360);v1(i)=quadv(fun1,0,y2(i)-2.05*tan(2*pi*4.1/360); v2(i)=quadv(fun2,y2(i)-2.05*tan(2*pi*4.1/360),0.4*tan(2*pi*4.1/360)+y2(i);endv3=v1+v2for i=1:3 y2(i)=1.18+0.01*(i-1); funa=(x)4.9*0.89*sqrt(1-(x-0.6).2/0.36);funb=(y)2*0.89*sqrt(1-(y-0.6)2/0.36)*(2.45-(y-y2(i)+2.05*tan(2*pi*4.1/360)/tan(2*pi*4.1/360);V1(i)=quadv(funa,0,y2(i)-2.05*tan(2*pi*4.1/360); V2(i)=quadv(funb,y2(i)-2.05*tan(2*pi*4.1/360),1.2);endV3=V1+V2模型三部分程序：clearclcsyms a y2=2.6;R=1.625;D=3;h=1;y1=2*tan(a)+y2;H=4*sin(a)*(sqrt(1.625)2-(y1-1.5)2)-0.625)/3*pi+y1;%if H=3 %H=3;%endt=2*acos(1-2*H/D);V1=0.5*t*(D/2)2*(h-R)-2/3*(R-h)3+2/3*R3)+sin(t/2)*cos(t/2)*(-(D/2)2*h+1/12*(D/2)4/R+1/120*(D/2)6/R3)+sin(t/2)*(cos(t/2)3*(1/6*(D/2)4/R+1/90*(D/2)6/R3)+1/45*(D/2)6/R3*sin(t/2)*(cos(t/2)5;y3=-6*tan(a)+y2;H=-4*sin(a)*(sqrt(1.625)2-(y3-1.5)2)-0.625)/3*pi+y3;t=2*acos(1-2*H/D);V2=0.5*t*(D/2)2*(h-R)-2/3*(R-h)3+2/3*R3)+sin(t/2)*cos(t/2)*(-(D/2)2*h+1/12*(D/2)4/R+1/120*(D/2)6/R3)+sin(t/2)*(cos(t/2)3*(1/6*(D/2)4/R+1/90*(D/2)6/R3)+1/45*(D/2)6/R3*sin(t/2)*(cos(t/2)5;syms yy1=2*tan(a)+y2;y3=y2-6*tan(a);fun1=2*sqrt(9-(3-y)2)*