运动界面追踪问题论文.docx

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

2、名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 参赛队员 (打印并签名) ：1. 孙元 2. 于冰 3. 陈魁东 指导教师或指导教师组负责人 (打印并签名)： 日期： 2010 年 7 月 16日赛区评阅编号（由赛区组委会评阅前进行编号）：2009高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：41 运动界面追踪问题孙 元于 冰陈魁东摘要本题专业性很强，运用到很多的专业理论公式，我们结合专业知识

3、和生活实际解决两个问题。对于问题一，我们针对管线的两种状态采用微元法进行分析。对于落水过程，我们建立了微分方程模型，模拟出水流下降过程，计算水流断面的速度和下降时间。对涨水过程，我们构建灌水模型，并分别讨论了前一峰点比后面的高和前一峰点比后面的低这两类情形，利用matlab计算出水流到达每一谷点和峰点的时间。计算结果为水流到达终点需要12.28天。对于问题二，我们在对有压管道非恒定流分析的基础上建立了水流的运动方程和连续方程，并利用特征线的计算方法将偏微分方程转化为常微分方程，然后对特征线方程组进行差分变换，配合初边值条件，编程求解输水管线全程的水头压力分布和从一种恒定流到另一种恒定流的时间。

4、由于在额定功率下最大输水流量为5.2m3/s；所以我们选取第一阶段和第二阶段流量的增量相等的方案；即第一阶段流量从0.6变为2.9；第二阶段从2.9升至5.2；这样我们得出两阶段的水压分布图以及总的时间关键词： 局部分析法 有压管壁模型 特征线法 初边值条件 问题重述某输水管线是利用114米地面落差有压、重力流输水工程。管径2.2米，管线长度176公里，地面高程变化很大。当输水管线建成后，首次通水时，为了防止水流速波动产生水击破坏管线，只能以每秒0.6立方米的流速由某水库向管线内灌水。当被供水城市水厂出流稳定后，再逐步提高入流速度达到每天45 万吨的额定工作状态。我们要解决两个问题1第一阶段：

5、假设以每秒Q（=0.6）立方米入流速度向管线内灌水，模拟出水头前沿面的运动过程，计算出水头前沿面顺序到达管线高程谷点和峰点的时间。2第二阶段：再分两阶段提高入流速度达到每天45万吨的额定工作状态，请设计通水方案，计算出由一种恒定流达到另一种稳定流后的管线水头压力分布和需要的时间。问题分析1. 问题一，我们针对管线的两种状态采用微元法进行分析。对于落水过程，我们建立了微分方程模型，模拟出水流下降过程，计算水流断面的速度和下降时间。对涨水过程，我们构建灌水模型，并分别讨论了前一峰点比后面的高和前一峰点比后面的低这两类情形，利用matlab计算出水流到达每一谷点和峰点的时间。 2. 问题二，我们在对

6、有压管道非恒定流分析的基础上建立了水流的运动方程和连续方程，并利用特征线的计算方法将偏微分方程转化为常微分方程，然后对特征线方程组进行差分变换，配合初边值条件，编程求解输水管线全程的水头压力分布和从一种恒定流到另一种恒定流的时间。由于在额定功率下最大输水流量为5.2m3/s；所以我们选取第一阶段和第二阶段流量的增量相等的方案；即第一阶段流量从0.6变为2.9；第二阶段从2.9升至5.2；这样我们得出两阶段的水压分布图以及总的时间术语和符号说明由于这道题是水力学问题，很多术语专业性较强，因此我们在此列出所用术语，有利于非专业人士阅读。符号物理量符号解释单位Q水流量单位时间内通过某一个过水断面的液

7、体的体积m3/sv断面平均流速单位时间内断面前进的距离m/sA过流断面面积过水段面是与所有流线正交的横断面m2c管壁湿周长度水流断面上流体与固体壁面接触的周界线mR水力半径过水断面面积与湿周之比mhf沿程水头损失单位重量的水体流动时由于边壁阻力在流程中所引起的水头损失mhj局部水头损失由局部边界急剧改变导致水流结构流速分布改变并产生漩涡区引起的水头损失mhw总水头损失沿程水头损失和局部水头损失之和mi水力坡度实际液体沿元流单位流程上的水头损失，即总水头线坡度1Z水流高程水流沿铅垂线方向到某水准面的距离me压力水头管壁对水的压强产生的水头mH管道水头压力即测压管水头，等于水流高程与压力水头之和m

8、根据以上定义结合物理规律，有以下公式成立： （1-1） （1-2） （1-3） （1-4） （1-5）注：P压强(帕)，水的密度(kg/m3)，g为重力加速度(m/s2)；符号解释同上。问题一在问题一中出水头前沿面的运动不外乎从高程下降、从低程上升和同高度向前传播三种状态。从所给数据可以发现绝大部分管线处于前两种状态，几乎不存在同高度管线，所以我们重点分析水流下降和上升两种状态。管线的总体高程变化图：模型假设：水流刚刚到达谷点时，谷点与相邻的峰点之间管线充水量很少，再后续计算时可以忽略。水流下降过程模型在最初充水阶段，由于水并未充满水管，因此水流下降时为无压流，我们取一小段管道中的水流结合动能

9、定理利用微元法分析： （1-6） （1-7）忽略项后化简得： （1-8）根据达西魏斯巴赫（DarcyWeisbach）公式： （1-9）式中，为无量纲待定系数，习惯上称为沿程阻力系数，根据已经定义的水力坡度表达式和谢才公式得： （1-10）再结合巴甫洛夫斯基公式 （1-11）式中，R为水力半径，y表达式： （1-12）本题中由于；所以近似计算时，y可取 （1-13）最后再讨论水力半径R的表达式：如上图，由几何关系可得：过水断面面积 （1-14） 湿周 （1-15）水力半径 （1-16）我们联立（1-1）（1-14）（1-15）（1-16）理论上可以整理得水力半径R和断面平均流速v的关系式，但实

10、际计算时由于R、V的关系式过于复杂，不易求解，我们利用matlab最小二乘法进行多项式拟合，近似求解，并得出表达式： （1-17）平均误差平方和为0.0184。拟合曲线与实际曲线对比图如下：将式(1-7)(1-10)(1-13)代人(1-8)式得： (1-18)其中R和V的关系式为式（1-17）由于利用matlab求解ODE方程。具体求解过程见附录程序。水流上升过程模型水流上升的过程中输水管线必然是充满水的，管壁对水流速的影响就不能忽略了，由于此时按照上述模型分析情况比较复杂，我们用总流量、流速以及时间之间的关系来求解水流上升所需的时间。问题一中流速Q一定，只需确定流量即可，根据连通器原理可知

11、该峰点与上一个峰点的绝对高度差对充水的总流量是有影响的，我们需要分别考虑两种情况：1）前一个峰点高于该峰点此时我们可得断水面上升到峰点时注入的水量体积为峰点高度以下的管线体积，如图所示，我们也可以认为峰点与相邻谷点的时间差为充满该部分水量的所需时间。体积=管线截面积管线长度时间=体积/流量Q2）前一个峰点低于该峰点此时由于断水面到达前一峰点时，有一部分管线已经充满水了，如图中横线所示。继续注水，直到水流充满峰点高度以下部分的管线，如图中涂黑部分以及竖线部分。由于节点与节点之间的距离较远，竖线部分水量较之涂黑部分微乎其微，可以忽略，所以我们可以近似认为峰点与相邻谷点的时间差为充满涂黑部分水量的时

12、间。体积=管线截面积管线长度时间=体积/流量Q模型的求解：我们通过编程利用迭代思想分段累加积分求解水断面到达各点的时间。假设入水口的初速度为V0=0.2，可以证明V0的大小对输水时间的影响不大。程序流程图：运行结果：峰点与到达峰点时间表：峰点时间80.02 420.35 520.47 580.72 660.99 701.15 721.47 761.53 791.59 831.76 951.94 1051.99 1102.10 1122.22 1162.37 1212.48 1252.63 1372.76 1432.85 1472.92 1493.01 1523.25 峰点时间1563.40 1

13、603.47 1633.53 1693.59 1823.76 1843.84 1863.90 1974.21 2004.23 2094.35 2144.45 2194.57 2254.65 2354.91 2404.99 2455.03 2505.10 2535.14 2555.17 2605.28 2625.37 2645.40 峰点时间2735.47 2785.59 2925.67 2955.75 3065.86 3126.04 3156.07 3226.14 3256.18 3356.23 3376.32 3426.36 3586.45 3686.75 3726.82 3786.89 3

14、886.98 4017.10 4087.31 4147.86 4178.45 4248.65 峰点时间4278.87 4349.00 4519.28 4689.78 47710.01 48210.24 49410.58 50010.71 51211.05 51511.14 52711.19 53911.32 54911.39 55511.45 56211.49 57311.55 58811.72 59411.91 59612.06 59912.10 60112.19 时间单位：天如下图：谷点与到达谷点时间表：谷点时间110.16 430.39 530.66 600.96 681.15 711.

15、46 731.52 771.57 801.75 921.94 961.96 1062.05 1112.22 1152.37 1182.46 1222.63 1322.74 1382.82 1442.91 1483.00 1513.25 1553.39 谷点时间1573.47 1613.52 1653.57 1743.74 1833.84 1853.90 1944.20 1984.23 2034.34 2104.43 2164.53 2204.61 2344.90 2374.98 2415.02 2495.10 2515.14 2545.17 2595.25 2615.37 2635.39 26

16、75.46 谷点时间2765.58 2825.66 2935.74 3025.84 3096.03 3136.06 3196.14 3236.17 3266.22 3366.32 3386.36 3456.42 3666.75 3696.82 3756.89 3836.98 3907.05 4037.23 4117.74 4158.36 4188.61 4258.78 谷点时间4298.96 4469.21 4629.77 47510.01 47810.14 49210.56 49710.69 51011.04 51411.14 51811.17 53311.32 54311.39 55311

17、.45 55911.48 56711.54 57711.69 59311.91 59512.06 59812.10 60012.17 60212.28 时间单位：天如下图：结果可知：水流传到终点需要12.28天。解决问题二有压管道非均流模型我们可以知道，在管道中的水流恒定的状态下，管道水流的截面积为一个圆形管道的截面，即水是充满整个管道的。可以从下图得出，水力学中，这种状态属于有压管道的状态。而管道中的水流由一种恒定流过渡到另一种恒定流的过程是一个非恒流的过程，因此我们使用有压管道的非恒流方程来求解。我们作如下假设：（1）压力管道的管壁被看作是弹性的。由于弹性水击更能反映水击的实际情况，因此我

18、们忽略刚性水击而只考虑弹性水击。（2）瞬变流的基本微分方程有运动方程和连续方程。运动方程从管道水体中选取控制体，应用牛顿第二定律可以推导出有压管道非恒流的运动方程。如下图：有压管道运动方程控制体在管道液体中选取长度为dx的微小控制体,x轴取与恒定流时的水流一致的方向，管轴线与水平线的夹角为a，则作用于微小控制提上的力如下图所示，上下游断面的水压力为，控制体周界面上的阻力为，侧水压力，以及重力mg。若上游面的密度为，过水断面的面积为A，湿周为，压强为P，则下游断面相应各量分别为，。则作用于微小控制体上的外力在x轴上的分力为【2】：(1) 上下游断面的水压力之差。 （2-1）(2) 控制体周界面上

19、的阻力（设控制体周边平均阻力为t ） （2-2）其中q 为控制体侧壁与管轴线的夹角，一般是很小的，可取.(3) 侧面水压力4 P 沿x轴的分量 （2-3）(4) 重力分量 （2-4）设控制体沿x轴方向的速度为v，则有压管道水体的加速度为 (2-5) 由牛顿第二定律可得，作用与x轴方向的所有外力的合力等于控制体的质量与沿x轴方向的加速的乘积，即： （2-6）因为v是时间t和坐标x的函数，所以可以得到： （2-7）把式(2-1)、(2-2)、(2-3)、(2-4)、(2-7)代入式(2-6)，并取整理可得： （2-8）因为测压管水头线，而控制体周边平均阻力可由达西魏斯巴赫公式表示，其中 f为恒定流

20、时的沿程阻尼系数，可将式（2-8）化为 （2-9）因为湿周，其中R为水力半径，将其代入式(2-9)可得： （2-10）此即有压管道非均流的运动方程。连续方程利用质量守恒原理可以直接推导出有压管流的非恒定流连续方程，在管路中选取两个非常接近的横截面作为控制体，两截面的间距为dx，如下图：我们给出最终推导结果： （2-11）此即有压管道非均流过渡过程连续方程。联立方程(2-10)，(2-11)，用特征线法求解，得到常微分方程组： （2-12） （2-13）在给出该方程的初始条件和边界条件后，即可通过迭代法求解第一个问题。我们的输水方案采取分两阶段输水，即先把流量从0.6m3 / s提升到一半即2.

21、9m3 / s，然后再提升另外一半至5.2m3 / s。水头压力分布的计算如下：首先我们假设有一种恒定流过渡到另一恒定流的过程中水头损失时局部水头损失可以忽略不计。沿程水头压力分布计算如下： （2-14）其中为入水口的水头压力，为管道上任一点的水头压力，i为水力坡度，x为管道上任一点距入水口的距离。代入=150，D =2.2 化解并整理可以得到： （2-15）再分别计算Q=2.9和Q=5.2时的压力分布：Q=2.9时：H0= 325.6此时式（2-15）化简为： （2-16）Q=5.2时：H0= 381此时式（2-15）化简为： （2-17）模型的评价模型优点：该问题是专业性很强的实际应用问题

22、，构建模型时我们充分考虑了输水的实际情况，由于输水管线较长，节点间距都是以千米度量，所以模型中计算体积时直接利用管线截面积与管线长乘积来表示是合理的，并且这样忽略弯头处的体积也方便计算。第一问，下降阶段利用微分方程精确的描述了水流的运动过程，但是在建立方程的过程中遇到一些项比较复杂无法计算，我们采用matlab 的拟合工具箱，化简方程，从而很好的解决了问题。全文中有丰富的理论基础，并配合大量的图像，从直观上和理论上给出了较为正确的结果。模型不足：模型中实际考虑出的物理因素太少，与实际输水问题还有一定差距，比如说我们没有考虑水库出水口的压力变化对出水口水压的影响；还忽略了很多实际的阻力因素，考虑

23、的仅仅是较为理想的输水管道问题，在一些专业的理论运用上，该模型也显得较为理想化。模型中还采用了大量的近似计算，在输水管线很长时这种近似还是可以接受的。参考文献1 吴持恭.水力学.高等教育出版社.20082 张德丰.MATLAB程序设计与典型应用.北京：电子工业出版社，2009.63 曲世琳.长距离输水管线的非恒定流动分析.【J】中国给水排水，2005.124 杨敏等.有压管道充水过程数值模拟.【J】水利学报，2007年02期5 穆祥鹏.长距离输水系统的过渡过程数值计算及水力特性研究.【J】天津大学建筑工程学院，2004年12月6 吴家鸣等，流体力学，【M】机械工业出版社，2005年11月附录第

24、一问程序：data1;data;p=data1(:,3);k=1;l=1;m=1;n=1;downnode=;upnode=;fengnode=;gunode=;downnode(l)=1;gutime=;fengtime=;l=l+1;s=zeros(1,602); %每次灌水的长度t=zeros(1,602); %流经两个节点的时间e=zeros(1,602); %各节点的局部阻力系数v=zeros(1,602); %水进入各节点的速度v1=zeros(1,602); %水流出各节点的速度tt=zeros(1,602); %水流到各节点所需时间w=zeros(1,5000); %水在两个节

25、点之间的速度node=;c=0.003549 0.1 0.003549 0.2 0.0894 0.003549 0.000588 1.5 0.00833 0.5 0.06 0.1 1.0 0.1;node(1)=1;global sinafor i=2:602; %确定各点属性 if (p(i)=p(i-1)&(p(i)p(i-1)&(p(i)=p(i+1) gunode(n)=i; n=n+1; end if p(i)p(i-1) upnode(k)=i; %上升节点 node(i)=2; for j=i:-1:1 %求与该节点等高的前一个点到该节点的距离 h=abs(data1(j,3)-

26、data1(j-1,3); L=data1(j,2)-data1(j-1,2); if data1(j-1,3)data1(i,3) h=abs(data1(j-2,3)-data1(j-1,3); s(k)=sqrt(h2+L2)+s(k); %距离累和 renode(k)=j; %等高节点（近似） break; end s(k)=sqrt(h2+L2)+s(k); end k=k+1; endendk=1;l=1;for i=1:602 if(i=1) h=data(1,3)-data1(1,3); L=data(1,2)-data1(1,2); sina=h/sqrt(h2+L2); w

27、(1)=1; for j=2:data1(i,2)-data(i,2) %划分两个节点间的距离 R=(34.4580*w(j-1)2+61.2317*w(j-1)-6.3853)/(358.934*w(j-1)2-29.8083*w(j-1)+2.6058)1.313; w(j)=w(j-1)+(sina-0.0012*w(j-1)2/R)/(0.102*w(j-1); %w=w+dw 近似微分 t(i)=t(i)+1/w(j); end tt(i)=t(i); v1(i)=w(j); v(i)=v1(i)*sqrt(1-e(i); %局部阻力影响 else for j=1:14 if(data(i+1,5)=j) e(i)=c(j); if(e=1) e=0.95; end