附录一：论文格式规范及论文模板.doc

2012年武汉科技大学数学建模竞赛
论文格式规范
l 论文（答卷）用白色A4纸单面打印，上下左右各留出至少2.5厘米的页边距。
l 论文第一页为控制页（附录二），论文题目和摘要写在论文第二页上。从第三页开始是论文正文。论文从第二页开始编写页码，页码必须位于每页页脚中部，用阿拉伯数字从“1”开始连续编号。注意，论文一律要求从左面装订。
l 论文不能有页眉，论文中不能有任何可能显示答题人身份的标志。
l 论文题目用三号黑体字、一级标题用四号黑体字，并居中；二级、三级标题用小四号黑体字，左端对齐（不居中）。论文中其他汉字一律采用小四号宋体字，行距用单倍行距，打印时应尽量避免彩色打印。
l 提请大家注意：摘要在整篇论文评阅中占有重要权重，请认真书写摘要（注意篇幅不能超过一页）。阅卷组评阅时将首先根据摘要和论文整体结构及概貌对论文优劣进行初步筛选。
l 引用别人的成果或其他公开的资料（包括网上查到的资料）必须按照规定的参考文献的表述方式在正文引用处和参考文献中均明确列出。正文引用处用方括号标示参考文献的编号，如[1][3]等；引用书籍还必须指出页码。参考文献按正文中的引用次序列出，其中：
Ø 书籍的表述方式为
[编号] 作者，书名，出版地：出版社，出版年。
Ø 参考文献中期刊杂志论文的表述方式为
[编号] 作者，论文名，杂志名，卷期号：起止页码，出版年。
Ø 参考文献中网上资源的表述方式为
[编号] 作者，资源标题，网址，访问时间（年月日）

论文模板（武科大全国数学建模竞赛2010年A题论文）
本文只是提供一种竞赛论文写作的大体格式，具体细节内容可视情况而定，并不一定要完完全全按照该论文的模式来写论文。大家关注的格式应该是以下这些：摘要 关键字 问题重述与分析 对题目某些条件的合理假设及变量说明 模型建立及求解 模型的评价与推广 参考文献 附录（主要是一些程序）

储油罐的变位识别与罐容表标定
摘要
本文研究了两种形状的储油罐在罐体发生变位后，变位角度参数的识别和罐容表的标定问题建立了比较精确的数学模型。
对于问题一，针对小椭圆形储油罐无变位的情况，通过微元分析法得到储油体积与油位高度的函数关系表达式（见正文3），并作出相应的曲线图像；通过与实际检测数据的比较和误差分析，可知无变位情况下所建立函数模型具有很高的精度。
在此基础上，针对纵向倾斜角度为的变位情况，我们也建立了储油体积和油位高度的函数模型：根据油位高度的不同，确定边界条件，分为三种情况进行讨论，得到了储油体积和油位高度见的函数关系，该函数是一分段函数（见附录2、附录3）。作出纵向倾角时的曲线图像，与实际数据的散点图比较检验，通过误差分析验证所得到模型的准确性；并计算出罐体变位后，油位高度间隔为1cm的罐容表标定值（见表一）。
对于问题二中的实际储油罐，我们按照以下步骤建立罐内储油量与油位高度及变位参数之间的函数关系：
首先分析仅有纵向偏转时罐内储油量与油位高度的函数关系。当纵向倾斜角度为时，根据罐内油位高度的不同，应分为三种情况进行讨论，通过几何分析的方法可得到该函数模型为（见正文）
再分析仅有横向偏转时储油量与油位高度的函数模型。对于横向偏转角度为的情形，通过几何关系，易得出实际油位高度与测量油高的关系式：，(为测量的油高，为圆柱体半径)。将代入无偏转时储油量的计算公式中，则可建立横向偏转时储油量的数学模型，（见正文）
在第三步中，我们假设储油罐先横向偏转度，再纵向偏转度，则在以上两步的分析基础上，可求出罐内储油量与油位高度及变位参数间的函数模型，。
根据建立的数学模型，我们采用穷举搜索法，计算理论值与实际检测值最小时的偏转角度，从而确定油罐的变位参数，，；然后根据公式求出了变位后油位高度间隔为10cm的罐容表标定值。我们利用所给的实际检测数据对所建立模型进行了比较和误差分析，发现理论数据与检测数据间的误差不大于3%，表明所建立模型具有较高的准确性。

关键字 微元分析法 曲线拟合 误差限 分段函数 比较检验

一、问题重述
已知加油站都有若干地下储油罐，且有与之配套的“油位计量管理系统”，通过预先标定的罐容表进行实时计算，以得到罐内油位高度和储油量的变化情况。
然而储油罐在使用一段时间后，由于地基变形使罐体位置发生纵向倾斜和横向偏转等变化，从而导致罐容表发生改变。因此需要定期对罐容表进行重新标定。题中给出了一系列示意图，并得出了一些实验数据和实际检测数据。
本题需要用数学建模方法研究解决储油罐的变位识别与罐容表标定问题。
（1）为了掌握罐体变位后对罐容表的影响，利用小椭圆型储油罐，分别对罐体无变位和倾斜角为的纵向变位两种情况做了实验，得出实验数据。请建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为的罐容表标定值。
（2）对于实际储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度和横向偏转角度）之间的一般关系。利用罐体变位后在进/出油过程中的实际检测数据，根据所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为的罐容表标定值。进一步利用实际检测数据来分析检验模型的正确性与方法的可靠性。

二、对问题的初步分析
要解决储油罐变位识别与罐容表标定问题，那么就需确定罐内油位高度与储油量的对应关系。第一问由于只有纵向倾斜，因此需要通过积分来分别求无变位和纵向变位后的油位高度与储油量的关系（，），最后通过对模型的求解来标定罐容表。而实际上，储油罐并不是如第一问那样简单的椭圆柱体，而且对于第二问不仅要考虑纵向倾斜，同时还要横向倾斜，因此我们要求出体积与油位高度的关系，必须要先确定。再对模型求解，标定罐容表。

三、问题假设
（1）假设储油罐是理想的，规则的罐体。不考虑物理形变，温度或其他因素导致的罐体形状不规则或储油体积的变化。
（2）假设油位探针、注油口、出油管等伸进油罐内的部分忽略不计。
（3）对于油罐左右倾斜求解方法相同，只因为油浮子位置导致计算结果有所改变，我们假设油罐按图中方位倾斜。
（4）假设图中标注的度量均不考虑油罐的厚度。

四、模型的建立及求解
符号规定：
：纵向倾斜角度 ：横向倾斜角度
：罐体左端液面高 ：椭圆长半轴
：椭圆短半轴 ：储油量
：对微元积分的微元高度 ：储油横截面面积

模型一：
研究罐体变位后对罐容表的影响，因此我们考虑先未变位时理论油位高度与储油量之间的关系，建立相应的数学模型。

设横截面椭圆的方程为：
（1） 先讨论无变位模型：
椭圆弓形的高为，图中带阴影部分为储油横截面，先用定积分求储油体积。
设椭圆弓形的面积为，则：

油罐的长为,储油的体积为，可得：

已求出储油量和油位高度的关系，用MATLAB,对的关系式作图，与实验数据作图进行比较（程序见附录1），得下图。

无变位进油量的数据图
红色的为曲线图
蓝色的为实际数据散点图

无变位出油量的数据图
红色的为曲线图
蓝色的为实际数据散点图

由图发现，未变位情况下，理论模型与实验数据几乎吻合，因此模型可用。

现在讨论变位后的模型。
因为当油罐纵向偏转角后，需分三种情况讨论。
1） 当油面到底面的投影时，如图：

对每一个椭圆部分面积微元在上进行积分，求体积；

∴

进行变量代换 则 ，

积分变为：

再用MATLAB求积分，由于积分结果较复杂，所以将计算结果及方法置于附录2。

2） 当油面将底面全都覆盖，又不到达顶面时（如图），及（单位为）的时候：

椭圆的面积

则是
令

积分变为：

再用MATLAB求积分（计算结果及方法同上见附录3）。

3） 当油面开始慢慢覆盖顶面，即倾斜时高端液面高恒定为1.2时（如图）：

同1) 设 进行变量代换得积分

MATLAB求积分同1） （计算结果及方法见附录2）
∴

而实验数据是从0.159开始进油，而从第二个数据开始，就完全进入油面将底面全都覆盖的情况，即讨论的第二种情况，因此，我们将建立的第二种情况的模型，与实际数据作图。

为了使作图简单，我们将油位高度全转化成用倾斜时高端液面高来表示

油位高度 （）用MATLAB得到两种曲线的图如下（作图方法见附录3）（横坐标为倾斜时高端液面油位高度）：

有变位进油量的数据图：红色的为曲线图，蓝色的为实际数据散点图

有变位出油量的数据图：红色的为曲线图，蓝色的为实际数据散点图

分析图像发现，模型做出来的曲线图与实际数据散点图几乎完全拟合。同时我们求临界点的数据，对于1），2）模型分别用求临界点时的体积，
得 ，

对于2）中曲线数据与散点数据用MATLAB（算法见附录4）求得
误差的最大值
误差的最小值
对于2）3）的临界点代入模型求得临界点

综上所述 误差范围都在最大误差0.0519内，因此，通过数据检验，模型正确。
最后，我们将三种情况的数学模型用MATLAB作图，将三段连接，得下图：

即为倾斜变位后油标高度与储油量的关系曲线。

现根据所建模型 用MATLAB(计算方法见附录5) 得出从油位高度为0开始 每隔1的罐体变位后罐容表的标定值，如下表：

表一

模型二：
对于实际储油罐，要建立罐体变位后的数学模型，我们先求
1）无变位情况即时实际储油罐体积与油位高度的关系

无变位的储油体积
油罐体积即求图中阴影部分的体积；
为了计算简单，先计算如图3所示灰色部分的体积。

体积计算的模型分两部分：
（1计算中间圆柱体部分储油的体积。
这部分相当于圆柱体的一部分。

只需要计算截面面积S。利用微积分，
有

则只需要计算截面,利用微积分：

（2）计算两端球缺部分所剩余油料的体积 2V
由于左右是对称的，所以我们只需要计算其中一个。

先计算水平截面的面积，然后在垂直方向上对积分，然后在垂直方向上对
的积分。

下面求：S

如图6，球的水平截面是一个圆

水平截面的半径

如图，就是阴影部分的面积
S=

=

＋

于是得到如图所示部分的体积计算公式

下面根据图3中阴影部分的体积来导出储油罐中体积和高度的关系式

对于题中的问题，由简单的几何关系可得r=1.625m;
根据推导的上式，可求得
最后代入实际问题的数据，可求得表达式（malab编程指令及计算结果见附录6）

2）变位后的实际储油罐体积与油位高度的关系
同样分三段考虑
（1）对于油面没有覆盖低端的情况如下图：
图①
设倾角为，油浮高为时，则
设圆柱体内油的体积为，设球冠ABC的体积为，球冠内BCD部分的体积为 则此时油的体积
设球冠体的半径为，圆柱的半径为，则
高为时，平行于圆柱底线的平面截球冠所得的圆的半径
所截弦长
如图：
图②
在上图中，用垂直于面的平面截该立体，则所得的平面如下图所示：
图③
面ABC和面的夹角为　∴　高
∴
为了求图②中立体体积　沿方向进行积分，建立坐标系如下：

则的坐标为
积分时，而

∴　　而　　　
∴
∴
令得　∴
∴原始的不定积分为
再令
∴上述积分为

∴
如果有横向偏转，则此时油的实际高度
那么横向偏转角度为，纵向倾斜角度为时，油的体积

（具体MATLAB计算程序见附录8）
第（2）种情况

阴影部分体积为须求的储油量，利用割补法得
又进行割补后，为无变位规则体型下的体积。所以用上述对无变位体积求得的公式，可得的表达式（见上对无变位的求解结果），又因为油浮高度
如果有横向偏转，则此时油的实际高度
将相关关系式代入表达式可得体积公式
第（3）种情况图形如下：

体积即为，这里为情况（1）种的体积。算法及结果见如上情况（1）
我们采用穷举搜索法，让、从到按步长0.01进行循环，即，取等于给出的实际测量油位高度值，计算，把所求得与给出的测量储油量数据进行比较求误差，计算理论值与实际检测值误差最小时的偏转角度，从而确定油罐的变位参数，（具体计算程序见附录8）。

4） 模型的分析与检验
现根据所建模型，我们画出实际数据的散点图及关系图如下（具体matlab编程指令及表达式见附录6）：

由图可看出，模型曲线和实际数据几乎完全重合，所以，模型可行性较好，且计算方法是可靠的。
因此用MATLAB(计算方法见附录7)得出从油位高度为0开始 每隔10的罐体变位后罐容表的标定值，如下表：

表二

五、模型的评价
1.模型的优缺点
（1）本模型数学推导严谨，理论性强。
（2）本模型充分利用MATLAB求积分、模型求解。作散点图，对MATLAB运用灵活，使
模型求解结果更准确，数据更可靠。
（3）模型建立中，算法、图形、数据各种表达方式充分利用，使模型更易理解，同时也很好的通过数据证明了模型的可行性。
（4）油罐变位后的第1）、3）段求出的模型，虽然与第二段临界处数据进行计算对比，在误差范围内，但没有实际数据做比较，因此对该段模型的正确性的评判不够全面。
（5）由于对模型进行求解时，结果比较复杂，而且难以给出简化模型，所以模型的实用性受到约束。
2．模型的改进
我们的模型数学推导严谨，完全用理论求解，最终也使得模型相对复杂，因此我们可以提出简化模型，利用多项式逼近，将复杂的积分结果转化成多项式，使模型的实用性更强。再考虑到测量数据导致的误差，为了使计算结果更准确，可以对模型一中的数据进行修正。另外，实际中大量的油料体积随着温度的改变有一定的变化。为了消除温度的影响，我们可以考虑油料的体积随温度变化的关系。可采用经验公式：，这样，我们可以将体积全部化为温度20度时的数据，然后进行计算。

六、参考文献
[1] 邓薇，MATLAB函数速查手册，北京，人民邮电出版社，2008。
[2] 李德宜、李明，数学建模，北京，科学出版社，2009。
[1] 孙宏达、关进波，用逼近法计算横截面为椭圆形（圆形）储油罐的储油体积，管道技术与设备，第3期，P29，2001。
[2] 潘孝光，倾斜卧式油罐容积测量与计算，油气储运，第6卷第6期，P48，1987。
[3] 田铁军，倾斜卧式罐椭直圆筒部分容积的近似计算，现代测量与实验室管理，2004年第1期，P1，2004。

附录1
无变位进油
%绘图的函数 m脚本文件
function huitu(x1,y1)
scatter(x1,y1,'b*')
hold on
x=0.15:0.01:1.19;
y=(pi/2+(x-0.6)./0.6.*sqrt(1-(x./0.6-1).^2)+asin(x./0.6-1))*0.6*0.89*2.45;
plot(x,y,'r-');
xlabel('油位高度h(米)');
ylabel('罐内总油量v(立方米)');
grid on;
box on;

无变位出油
%多项式拟合
p=-1:0.01:1;
q=pi/2+p.*sqrt(1-p.^2)+asin(p);
q1=pi/2+16/(1575*pi)*(615*p-80*p.^3-48*p.^5);
plot(p,q1,'r*')
hold on;
xlabel('x ');
ylabel('y');
plot(p,q,'k.');
grid on;
box on;

附录2
%求定积分
sym t
f='0.5*pi+t*sqrt(1-t^2)+asin(t)'
V=int(f,'t',-1,'h/b-1')
V =asin(h/b - 1.0)*(h/b - 1.0) + (1.0 - 1.0*(h/b - 1.0)^2)^(1/2) - 0.33333333333333333333333333333333*(1.0 - 1.0*(h/b - 1.0)^2)^(3/2) + (1.5707963267948966192313216916398*h)/b - 1.5707963267948966192313216916398

附录3
Malab求积分

syms a b h r ;
int(‘pi/2+x*sqrt(1-x^2)+asin(x)’,’x’,’h/b-1-2.45*tan(r)/b-1’,’h/b-1’)
V(x)= (a*b^2*((1 - (1/b*h - 1).^2).^(1/2) - (1 - (1/b*h - 1).^2).^(3/2)/3 - (1 - (49/20/b*tan(r) - 1/b*h + 1).^2).^(1/2) + (1 - (49/20/b*tan(r) - 1/b*h + 1).^2).^(3/2)/3 + asin(h/b - 1).*(h/b - 1) + (pi*(h/b - 1))/2 + (pi*((49*tan(r))/(20*b) - h/b + 1))/2 - asin((49*tan(r))./(20*b) - h/b + 1).*((49*tan(r))./(20*b) - h/b + 1)))/tan(r)
% 第二段曲线的拟合情况

a=0.89;
b0.6;
h=0.4:0.01:1.1;
r=pi/180*4.1;
v=(a*b^2*((1 - (1/b*h - 1).^2).^(1/2) - (1 - (1/b*h - 1).^2).^(3/2)/3 - (1 - (49/20/b*tan(r) - 1/b*h + 1).^2).^(1/2) + (1 - (49/20/b*tan(r) - 1/b*h + 1).^2).^(3/2)/3 + asin(h/b - 1).*(h/b - 1) + (pi*(h/b - 1))/2 + (pi*((49*tan(r))/(20*b) - h/b + 1))/2 - asin((49*tan(r))./(20*b) - h/b + 1).*((49*tan(r))./(20*b) - h/b + 1)))/tan(r);
plot(h,v,'r-');
hold on;
p=[0.439962356 0.452122356 0.467002356 0.479212356 0.492572356 0.506412356 0.518042356 0.531232356 0.543362356 0.555512356 0.567552356
0.580632356 0.593072356 0.605232356 0.617412356 0.628232356 0.640292356
0.652112356 0.664252356 0.674952356 0.687262356 0.698892356 0.709302356
0.721702356 0.733342356 0.745122356 0.756332356 0.768062356 0.779572356
0.790222356 0.802102356 0.814062356 0.824712356 0.836942356 0.849472356
0.861472356 0.873142356 0.884962356 0.896272356 0.908732356 0.921592356

0.933012356 0.946012356 0.958572356 0.970092356 0.983272356 0.996762356
1.008812356 1.021082356 1.035012356 1.047742356 1.062912356 1.064032356
];
q=[0.96286 1.01286 1.06286 1.11286 1.16286 1.21286 1.26286 1.31279 1.36279 1.41273 1.46273 1.51273 1.56273 1.61273 1.66273 1.71273 1.76273
1.81273 186273 1.91273 1.96273 2.01273 2.06273 2.11273 2.16273 2.21273
2.26273 2.31273 2.36273 2.41273 2.46273 2.51273 2.56273 2.61273 2.66273
2.71273 2.76273 2.81273 2.86273 2.91273 2.96273 3.01273 3.06273 3.11273
3.16273 3.21273 3.26273 3.31273 3.36273 3.41273 3.46273 3.51273 3.51474
];
scatter(p,q,'b*')
xlabel('油位高度h(米)');
ylabel('罐内总油量v(立方米)');
grid on;
box on;

附录4
求误差
function wucha(h,v)
a=0.89;
b=0.6;
r=pi/180*4.1;
v1=(a*b^2*((1 - (1/b*h - 1).^2).^(1/2) - (1 - (1/b*h - 1).^2).^(3/2)/3 - (1 - (49/20/b*tan(r) - 1/b*h + 1).^2).^(1/2) + (1 - (49/20/b*tan(r) - 1/b*h + 1).^2).^(3/2)/3 + asin(h/b - 1).*(h/b - 1) + (pi*(h/b - 1))/2 + (pi*((49*tan(r))/(20*b) - h/b + 1))/2 - asin((49*tan(r))./(20*b) - h/b + 1).*((49*tan(r))./(20*b) - h/b + 1)))/tan(r);
dv=v1-v;
m=dv./v1;
max(m)
min(m)
结果
ans=
0.0519 误差 最大值
ans =
0.0156 误差 最小值
附录5
%求罐容表的标定值
h1=15:1:117
h1=h1/100;
a=0.89;
b=0.6;
r=pi/180*4.1;
h=h1+0.4*tan(r);
v=(a*b^2*((1 - (1/b*h - 1).^2).^(1/2) - (1 - (1/b*h - 1).^2).^(3/2)/3 - (1 - (49/20/b*tan(r) - 1/b*h + 1).^2).^(1/2) + (1 - (49/20/b*tan(r) - 1/b*h + 1).^2).^(3/2)/3 + asin(h/b - 1).*(h/b - 1) + (pi*(h/b - 1))/2 + (pi*((49*tan(r))/(20*b) - h/b + 1))/2 - asin((49*tan(r))./(20*b) - h/b + 1).*((49*tan(r))./(20*b) - h/b + 1)))/tan(r);
v=v*1000

附录6
r=1.625;
R=1.5;
L=8;
v=pi*R^2*L+2/3*pi*(r-sqrt(r^2-R^2))*(2*r+sqrt(r^2-R^2))
求
syms a h r R L;
r=1.625;
R=1.5;
L=8;

v=pi*h*r^2-pi/3*h^3+h*L*sqrt(R^2-h^2)-4/3*h*sqrt(r^2-R^2)*sqrt(R^2-h^2)+(L*R^2-2/3*sqrt(r^2-R^2)*(2*r^2+R^2))*atan(h/sqrt(R^2-h^2))+4/3*r^3*atan(h*sqrt(r^2-R^2)/(r*sqrt(R^2-h^2)))+2/3*(h^3-3*h*r^2)*atan(sqrt(r^2-R^2)/sqrt(R^2-h^2));
结果如下：
v=(247*pi)/24 + (5707*atan((h - 3/2)./(9/4 - (h - 3/2).^2).^(1/2)))/384 + (2197*atan((5*(h - 3/2))./(13*(9/4 - (h - 3/2).^2)^(1/2))))/384 + atan(5./(8*(9/4 - (h - 3/2).^2).^(1/2))).*((2*(h - 3/2).^3)/3 - (169*h)/32 + 507/64) + (43*(9/4 - (h - 3/2).^2).^(1/2).*(h - 3/2))/6 + (169*pi*(h - 3/2))/64 - (pi*(h - 3/2).^3)/3;
绘图的编程指令：
function huitu1(x,y)
scatter(x,y,'g+');
hold on;
for h=0:0.05:3;
t=(247*pi)/24 + (5707*atan((h - 3/2)/(9/4 - (h - 3/2)^2)^(1/2)))/384 + (2197*atan((5*(h - 3/2))/(13*(9/4 - (h - 3/2)^2)^(1/2))))/384 + atan(5/(8*(9/4 - (h - 3/2)^2)^(1/2)))*((2*(h - 3/2)^3)/3 - (169*h)/32 + 507/64) + (43*(9/4 - (h - 3/2)^2)^(1/2)*(h - 3/2))/6 + (169*pi*(h - 3/2))/64 - (pi*(h - 3/2)^3)/3;
scatter(h,t,'r>');
hold on;
end;
xlabel('油位高度h(米)');
ylabel('罐内总油量v(立方米)');
grid on;
box on;

附录7
h=0:0.01:3;
f=(247*pi)/24 + (5707*atan((h - 3/2)./(9/4 - (h - 3/2).^2).^(1/2)))/384 + (2197*atan((5*(h - 3/2))./(13*(9/4 - (h - 3/2)^2).^(1/2))))/384 + atan(5/(8*(9/4 - (h - 3/2).^2).^(1/2)))*((2*(h - 3/2).^3)/3 - (169*h)/32 + 507/64) + (43*(9/4 - (h - 3/2)^2).^(1/2).*(h - 3/2))/6 + (169*pi*(h - 3/2))/64 - (pi*(h - 3/2).^3)/3;

附录8
h0=h0/1000; %h0为给定罐容表中的油位高度
v0=v0/1000; %v0为给定罐容表中的储油量
求误差最小时的偏转角度
% h0=h0(50:100);
% v0=v0(50:100);
r=1.5;
angA0=0;
angB0=0;
error0=10000;
for angB=0:0.1:5
u1=(angB/180)*pi;
h=r+(h0-r)*cos(u1);
for angA=0:0.1:5
u2=(angA/180)*pi;

for i=1:length(h)
v(i)=QingxieV(h(i),u2);
end

error=sum((v-v0').^2);
if error angB0=angB;
error0=error;
end

end

end
angA0
angB0
error0
%%计算倾斜后带球冠的油罐的装油体积，先算油比较低的时候
function vv=QingxieV(hh,angA)
%clear all;
syms x;
r=13/8;R=1.5;
h=hh;
ang0=angA;
u=(ang0/180)*pi;

%h=0.42;ang0=3;
%v1=YuanzhuV12(h,ang0);
%v1=YuanzhuV1(h,ang0); 求圆柱里面油的体积
if h<=0.42
v1=noChange2(h); %*********

h2=h+3*tan(u);
v2=qiuguanV(h2);
%v2=qiuguanV(h2);

%计算多加的那一部分油的体积
h3=h+3*tan(u);
%h3=h;
rj=sqrt(2*r*h3-(h3)^2);

x1=sqrt(2*r*h2-(h2)^2)-sqrt(2*R*h2-(h2)^2);
x2=rj;
%x2=r;

a=rj;

%t=rj-x1;
t=sqrt(rj^2+rj-sqrt(rj-x1));

pp1=- t^9/9 + ((3*a^2)/7 + (3*a)/7)*t^7 + (- (3*a^4)/5 - (6*a^3)/5 - (3*a^2)/5)*t^5 + (a^6/3 + a^5 + a^4 + a^3/3)*t^3;
v31=(tan(u)+tan(u)*sec(u))*pp1;

%p=rj-x2;
p=sqrt(rj^2+rj-sqrt(rj-x2));
pp2=- p^9/9 + ((3*a^2)/7 + (3*a)/7)*p^7 + (- (3*a^4)/5 - (6*a^3)/5 - (3*a^2)/5)*p^5 + (a^6/3 + a^5 + a^4 + a^3/3)*p^3;
v32=(tan(u)+tan(u)*sec(u))*pp2;

v3=v32-v31;
zongV=v1+v2-v3;
vv=zongV;
end

%%*******************
hL=h-3*tan(u);
hH=h+3*tan(u);
if h>0.42
myv1=zongV2(hL);
myv2=zongV2(hH);
myv=(myv1+myv2)/2;
vv=myv;
end
%myv

%%计算没有偏转的时候圆柱的装油体积
function VnoChange2(d0)
syms x y S V0
a=1.5; b=1.5; L=8;
d0=d0;
%d0=1.19349;
%d0=0.15902;
if d0<1.5
d=d0;
pp=0;
else
d=3-d0;
pp=1;
end

y0=b-d;
x0=(a/b)*sqrt(2*b*d-d^2);

y=(b/a)*sqrt(a^2-x^2);
S=2*int(y-b+d,0,x0);
V0=S*8;

if pp==0
V0=V0;
else
V0=pi*a*b*L-V0;
end

V=eval(V0);
end

25