数学建模校赛第三次作业

1、西南交通大学2014年大学生数学建模竞赛题目： B （填写A、B或C题）摘要供水管网是城市建设和发展的重要基础设施，是保障人民生活和发展经济建设的不可缺少的物质基础。在如今快速发展的社会中，广大民众对供水管网安全性的要求提高到了前所未有的高度。因此对于如何提高与监控城市供水管网安全就需要进行科学研究及分析。第一问主要涉及流体力学相关知识，第二问则是一个合理性检验问题，第三问需要动态监测水质污染情况，第四问则是对动态水质监测的规划问题。根据这些特征，我们对问题1用微元法建立单节点多管水流模型的方法解决；问题2我们采用配对样本的非参数检验法来研究样本观测量和指定的理论分布是否吻合；问题3我们借用医

2、学上药物动力学的二室模型进行求解；最后对于问题4我们在偏微分方程的基础上建立水管网水质动态模型对其进行模拟。问题1我们首先建立了双节点单管水流模型，接着我们运用流体动力学和牛顿第二定律，对水流进行微元受力分析，借助Matlab软件，对模型进行了求解，得出了节点流量与管道长度、地面标高间的关系；其次，在对双节点单管水流模型改进的基础上建立了单节点多管水流模型，即考虑一个节点的多进多出情况，最终得到在接近真实情形时的节点流量与管道长度、地面标高之间的关系。问题2是对问题1的补充验证，我们采用问题1中的单节点多管水流模型，对模型和数据采用配对样本的非参数检验法，运用SAS软件，对附件中所提供的数据进

3、行了筛选，去除异常数据,对残缺数据进行适当补充,并通过所有数据对理论结果进行了数据模拟，结果显示，概率，大于0.05置信区间的显著性水平，因此接受原假设，即实际值可以反映理论值，最后我们借助Epanet水利工程软件通过模拟运行进行了误差分析。问题3我们首先利用药物动力学，建立了房室模型，在对房室模型改进的基础上建立了二室模型。接着根据污染物质量守恒原理，建立了带参数的微分方程组。接下来借助Epanet水利工程软件，对苯的扩散污染进行了模拟，得到模拟的观测数据。然后我们运用Matlab中提供的最小二乘估计法，拟合求得微分方程组中的参数，从而求解出每个节点苯浓度方程为：。本文分别建立了从给水管网微

4、观模型到宏观模型，重点采用微分和偏微分数学建立方程，求解过程中用了一阶拟线性偏微分方程解法等多种微方程法法，用Matlab，Mathematica软件辅助求解方程，最终用Epanet水利软件模拟数据，对我们的估计值进行检验。关键词：单节点多水管模型，动态模型，非参数检验法，二室模型，Epanet，Matlab一、问题重述供水管网是城市建设和发展的重要基础设施，是保障人民生活和发展经济建设的不可缺少的物质基础。2014年4月10日兰州发生自来水苯含量超标事件。兰州市威立雅水务集团公司检测显示，4月10日17时出厂水苯含量、10日22时自流沟苯含量、11日2时自流沟苯含量均远超出国家限值的10微克

5、/升。2014年4月12日13:13原因已经查明：兰州自来水苯超标系兰州石化管道泄漏所致。直至2014年4月15日，兰州主城区的城关、七里河、安宁、西固四区已经全部解除了应急措施，全市自来水恢复正常供水。通过兰州自来水污染事件，广大民众对供水管网安全性的要求提高到了前所未有的高度。因此对于如何提高与监控城市供水管网安全就需要进行科学研究及分析。附件1第3节中给出了安徽某县城的供水管网的相关资料及数据，请结合上述数据资料回答以下问题：1、 构建数学模型，探讨节点流量与管道长度、地面标高间的关系；2、 利用附件1所提供的数据来验证上述模型的合理性；3、 若在节点17处出现了苯泄漏而导致该节点管道笨

6、含量超标的突发事件，请建立模型研究在不考虑取水时该物质在供水管网整体变动情况；4、 若要建立对该县城供水管网水质的监控点，以此能对城市管网的水质能进行实时监控，请确定监控点的设置数量及位置；5、 某日若发现了自来水中笨含量严重超标的情况，并收集到了第四问中你提出的各个监控点的数据资料，如何利用这些资料来估计笨的泄漏点。二、问题分析2.1 问题一的分析 问题一要求我们探讨节点流量与管道长度、地面标高间的关系，参考流体力学的相关内容，可以得到水头损失与流量、管段长度的函数关系，这里的水头损失即水流在水管中流动受到阻力损失的能量。首先建立双节点单管水流模型，即考虑水流在一根水管中从一个节点流向另一个

7、节点的情况，由此对供水管网中的水进行分析，用微元法取单位质量的小水快运用牛顿第二定律对其进行受力分析，建立微分方程组，运用Matlab求解，可以得到节点流量与地面标高、管道长度之间的函数关系。 其次，在第一个模型的基础上，考虑到实际情况，我们对它进行优化与改善，建立单节点多管水流模型，即考虑水流从多管道流向一个节点，再从该节点从多管道分流而出，因为水流从初始进入管道到达到平衡是一个复杂的过程，通过查阅相关专业文献，可得节点流量和沿线流量之间的相互转化关系，再将双节点单管水流模型中的结果看作是无线接近节点的沿线流量。由此，即可得在复杂供水管道模型中的节点流量与地面标高、管道长度之间的函数关系。2

8、.2 问题二的分析 问题二要求我们利用附件1所提供的数据来验证上述模型的合理性，即根据问题一中的机理关系，通过附件一给出的地面标高和管道长度的数据，求出各节点流量，再与原图中各节点流量的数据相比较，来检验问题一中模型的合理性。因为供水管道流量的影响因素有很多，而问题一中只是从一个方面来进行分析，因此理论值可能与实际值有较大出入。在这里，需要对数据进行拟合优度检验，我们采用两个配对样本的非参数检验法检验法来研究样本观测量和指定的理论分布是否吻合，即研究各点流量理论值之间的数值关系与各点流量实际值之间的数值关系是否相吻合。2.3 问题三的分析问题三要求我们对水管网各处污染浓度做出评测与判断。我们考

9、虑到污染源为17节点，根据附录1提供信息，我们将研究对象明确到1722号节点。要建立污染物扩散模型，首先我们要考虑到污染源在各节点之间的扩散和吸收过程，因考虑到污染物和各节点污染物扩散关系大致相同，所以关于污染物扩散浓度的计算可以使用医学上药物动力学的房室模型进行求解，本模型就采用的是二室模型对所求节点污染状况和污染源进行研究。最终我们就能对水管网中所有的节点浓度进行评测与判断。2.4 问题四的分析问题四需要我们选取一个最佳的观测点，然后对整个水管网的污染情况作出判断，我们选取污染变化速率为选取准则，当变化率最快时，能够最快的预测污染趋向并且提高测量的精确度。由于我们主要研究管道和水质变化，因

10、此我们选取管网水质动态模型对水管污染情况进行动态模拟。水流之间污染物浓度变化，我们采用一维质量守恒的偏微分方程描述，对动态模型，因为我们以流场中每一流体质点作为描述对象，所以我们采用拉格朗日时间驱动法对模型求解，其中运用一阶拟线性偏微分方程求解法。3、 模型假设l 本题中供水管道形状为空心圆柱体。l 苯泄漏过程只考虑苯泄漏后迅速与水混合扩散的情况，不考虑苯在水管壁上的流动过程。l 不考虑水流在节点处的湍流和涡旋。l 不考虑在供水系统中工业或居民集中供水的情况。l 苯不会与水中其他物质发生化学反应。l 不考虑其他因素对水管中水流的影响。4、 符号说明序号符号符号说明1水在水管内的流量2水在水管内

11、的流速，表示节点初速度3管道横截面积4摩擦阻力系数5节点的地面标高6水管长度7水管直径8污染源向其他节点的渗透度9其他节点对污染源的渗透度10污染物初始污染量11第i节点苯的污染浓度12苯浓度变化速率13比流量14介质密度（单位：）15集中流量（默认为0）16节点流量17供水系统总用水流量18初始污染物五、模型的建立与求解5.1 问题一的模型建立与求解5.1.1 问题一模型的建立 鉴于水在水管中的运动是复杂的，且受到多种因素的影响。因此，在问题一中，我们考虑的情况为双节点单管水流模型，具体阐述如下：如图5.1.1所示，节点A和节点B为单水管的两端节点，其中、之间的距离为，节点的标高为，节点的标

12、高为，水流流向为到。由于在供水管网节点处无外力作用，不妨假设节点的左侧和节点的右侧都不会出现分流的情况，因此对于节点和而言，节点流量即为单水管中、位置处的管道流量。因此，问题一在此模型下可以重新表述为：求点处的流量与、两点的标高、和、两点间的管道长度之间的关系，其中点的流量已知为。图5.1.1 为了求得三者之间的关系，可以由双节点单管水流模型下建立微分方程模型，具体叙述如下：根据流体动力学和输水管道中水流运动知识，水流在水管中满载且为不可压缩流，在水流从到源源不断地流动过程中，既不会挤压使体积变小，也不会出现断流情况，因此水管中水流运动情况都是相同的，所以我们考虑单位体积的水从到的运动过程，质

13、量为m。对所取水流微元进行受力分析：它受到自身重力在AB方向上的下滑分量、管壁对它的摩擦力和周围水对它的粘稠阻力。而水的流量和速度有如下换算关系： 其中为流量，为水流速度，为管道的横截面积。 因此，在水管的横截面积一定的情况下，水管中的水流流量可以由流速和横截面积的乘积来表示。 由受力分析写出它的运动方程为： （5.1.1） 其中，为管壁摩擦力和水粘稠阻力之和，为重力加速度，为水流微元的加速度。 由物理知识和查阅文献得到： 其中为水流的速度（单位：），为介质密度（单位：），为单位质量的水在管道内部的长度（单位：），为水力直径（即在满载状态下的水管直径，单位：)，为管壁对水流的摩擦阻力系数。 则

14、根据式5.1.1的运动方程写出微分方程模型为： 其中 。5.1.2 问题一模型的求解 运用Matlab求解微分方程组:将微分方程带入到Matlab中求解得到关于的解析解（） （5.1.2） （5.1.3）即得到关于时间的参数方程。 令，代入到式（5.1.3）中，反解出，即得到水流从点到点的时间，再将代入到式（5.1.2）中，即可得到速度与管道长度、两点的标高之间的关系。这时再另，即可得到节点流量与管道长度、地面标高间的关系，最终结果如下：其中，为初始条件。5.1.3模型一的完善与改进 在双节点单管水流微分方程模型中只考虑了单管且不分流的情况，但事实却往往并非如此简单。在实际问题当中，水管中的水

15、流节点通常意义上的概念即为“分流之点”，因此，在双节点单管水流微分方程模型的基础上，我们要对它进行进一步的改善，考虑到供水管道中复杂的汇总与分流情况，我们提出可以应用于实际问题的第二个模型：单节点多管微分方程模型。具体阐述如下： 图5.1.2 如图5.1.2所示，在图5.1.1的基础上，我们对节点加以考虑，在点的左右两侧分别加入了根入水管和根出水管，即水流分别从节点流向节点，再从节点分流到节点。这时，节点的流量就不单单是某一条入水管中的管道流量，而是这条入水管水流在点处汇总的结果。；而同样，对于出水管而言，它们的初始流量也并非的节点流量，而是这条出水管水流在点处的分流结果。因此，在双节点单管模

16、型中得出的点节点流量，相当于在单管多节点模型中每一条入水管道接近节点处的沿线流量（并未汇总）。1、输水管沿线流量向节点流量地转化： 通过查阅文献可知，对于一个复杂的输水管管道，每一个节点流量包括两部分：沿线流量化成节点流量和该节点的集中流量之和： 在不考虑居民或企业集中用水的情况下，集中流量。因此，的节点流量可以表示为： 其中表示管道接近的沿线流量，无限逼近之后即为双节点单管水流模型中点处的节点流量结果。因此，通过这个结果可以在输水管沿线流量已知的前提下，可以求出点的节点流量。2、节点流量向出水管沿线流量的转化： 在供水管道系统各节点流量已知的情况下，有如下关系： 其中表示比流量，表示供水系统

17、总用水流量，表示集中总供水量，表示总管道长度。而出水管的沿线流量可以表示为： 其中表示沿线流量，表示段的管道长度。通过上述两个方程可以求出出水管的初始沿线流量，再代入双节点单管水流微分方程模型可以求出的节点流量。因此，通过迭代的方法，在初始节点流量已知的前提下，可以依次求出每个节点的节点流量和每条管道的沿线流量。 由此，在一个比较接近实际的复杂网络供水模型中，我们可以得到每一处节点流量和地面标高以及管道长度之间的关系。5.2 问题二的模型建立与求解5.2.1数据的处理 首先，对图中数据进行处理和汇总；其次，将图中各节点地面标高和各管道长度带入到问题一的单节点多管微分方程模型中，根据从自来水厂流

18、出的初始流量，可以依次算出每一处的节点流量的理论值，并将它们绘到一张散点图中。 图5.2.15.2.2问题二模型的建立与求解因为问题二是让我们用附件中的数据来验证问题一中的模型，故附件中数据可以看做为已知的分布总体，而由问题一中模型得出来的结果看作是样本观测值。为了验证理论值与实际值，我们采用了两个配对样本的非参数检验法，利用SPSS软件进行数据分析，共分四个步骤进行。 （1）提出假设 这里我们定义原假设和被择假设分别为： 原假设：实际值可以反映理论值 备择假设：实际值不可以反映理论值 （2）计算统计量 我们通过SPSS软件给出各统计量的值，结果如图5.2.3所示： 描述统计量图5.2.3（3

19、） 根据给定的显著性水平的样本个数n，确定两个配对样本非参数检验的临界值（4） 若,则在的显著水平上，不能拒绝;否则，拒绝。5.2.3 误差分析 在实际操作中，因为水流在水管中的流动情况是相当复杂的，最终的结果是与实际值有较大误差的，所以在上述模型中我们仅仅考虑的是实际值与理论值的分布是否相同，而非数据值的完全相同。因为数据量少等原因的限制，我们无法利用其他方法对我们的模型数据进行误差分析。所以在这里，我们利用城市水管建设工业专用软件EPANET对水流流量进行了动态模拟，即将问题一单节点多管水流微分模型模型中求出的节点流量作代入到该软件中作为初始条件，而最终模拟的结果与我们用上述模型求出的结果

20、比较相符，因此我们认为问题一中的模型具有合理性。模拟结果如图5.2.4所示。图5.2.45.3问题三的模型的建立与求解5.3.1问题三模型的建立 问题三要求对泄漏物污染浓度进行监控。 对苯的泄漏点和检测点进行研究，泄漏点向下游进行了渗透和污染的，反过来下游对泄漏点也有小部分渗透，我们分别用，表示渗透度，其中表示泄漏量，表示浓度。 因此根据图中信息可写出两者苯浓度，满足的微分方程： （5.3.1）在逐渐向源注入污染物的时候，有：，。综合可得微分方程组模型如下： （5.3.2） 解一阶线性非其次微分方程，则它对应的其次方程通解为： 其中以由下列关系来确定：5.3.2问题三模型的求解5.3.2.1拟

21、合 我们运用Epanet软件，模拟19节点污染状况，得到模拟的观测数据运用Matlab中提供的最小二乘估计，得到5.3.2.2确定， 当时，因此： 由此确定： 将所有参数代入式，得到：通过对t的等步长取点得到散点图，再利用Matlab多项式拟合得到函数图像（各拟合参数详见附录四）：5.4问题四模型的建立与求解 由于问题所求为水质浓度的动态变化率，因此我们建立偏微分水质动态模型，采用拉格朗日时间驱动法解决水质动态模型。5.4.1问题四模型的建立 在水质动态模型中，在一个水力步长内,每一管道在某一时间段上的流速和流向已经确定,并假定每一管道的水流状态保持不变,水中物质的迁移和反应速率保持不变。在每

22、一连续水流片段中,每一水体的污染物首先发生反应,接着转移到下一单元水体中。当邻近水体为一节点时,节点处水流发生完全混合.管道中的反应完成以后,每一节点的物质混合浓度可以计算得出。因此，对于每一段水管中污染物浓度变化的偏微分方程为： 求解偏微分方程得到：，接着，对函数对时间微分，得到浓度变化率函数： 接下来对于浓度率变化函数，当取最大值时，说明在管道i中x处浓度变化率最大。即在x处，最易检测到浓度变化，即设为观测点。5.4.1问题四模型的求解 对于偏微分方程： 首先此方程的特征方程组： 接着得到首次积分： 然后我们发现该问题属于Cauchy问题，因此利用Cauchy问题求解方法得到： 从而求得：

23、 最后对每段求浓度变化速率函数：，令时，有：，此时，在i管道于处建立观测站，使得更容易观测得到数据。六、模型的评价与改进6.1模型的评价6.1.1模型优点总结（1）问题一中为了简化过程，建立了双节点单管水流微分方程模型，可以很好的求出沿线流量（可以看成节点流量）与地面标高、管道长度之间的关系；而为了比较切合实际，讨论复杂供水情况时，在上述模型基础上，又提出了单节点多管道水流微分方程模型，由简入深，由易入难，最终在比较符合真实情况的条件下，我们求出了节点流量和地面标高、管道长度之间的关系。（2）问题二中我们利用了SPASS和EPANET软件对两个较少的数据样本进行了非参数检验，最终结果证明了问题

24、一的合理性。（3）问题三中我们运用二室模型，较好的模拟了两节点间浓度的变化，并且能够从两节点推广到多节点浓度变化之间的关系。6.1.2模型缺点总结（1） 问题一中考虑的情形比较简单，没有考虑水流在汇流处的湍流和涡旋情况，并且将水流之间的粘稠阻力合并为摩擦阻力，有一定的误差。（2） 问题二中仅仅对23个节点流量数据进行了分析，显然数据量不够大，可能导致误差较大。（3） 题目张要求22节点最后无水池，但是根据情况，我们将第22个节点后添加了水池，以便浓度的扩散。6.2模型的改进（1） 考虑节点处的湍流和涡旋情况，并且考虑水流之间的粘稠阻力，根据流体动力学知识，可以得到复杂的水流模型。（2） 考虑污

25、染物渗透出水管道，污染物浓度将发生变化。7、 参考文献1董晓磊,沈致和.城市给水管网水质研究动态.合肥工业大学学报(自然科学版).(将于2006年第八期发表)2关于酒精在人体内扩散及残留时间问题的研究_陈宇3偏微分方程教程 （美）密苏里大学著4流体动力学 （美）W.F.休斯 J.A.布赖顿著5SAS统计分析实用宝典 清华大学出版社6流体动力学计算 中国原子能出版社7数学建模第四版 姜启源等 著附录：附录一：问题一微分方程的求解x,y=dsolve('k3*Dx=k1-k2*x2,Dy=x','x(0)=v0,y(0)=0')f=sym('-1/2/k2*

26、k3*log(tanh(k1*k2)(1/2)*(t+1/2*log(-(k1*k2)(1/2)+v0*k2)/(-(k1*k2)(1/2)+v0*k2)*k3/(k1*k2)(1/2)/k3)-1)-1/2/k2*k3*log(tanh(k1*k2)(1/2)*(t+1/2*log(-(k1*k2)(1/2)+v0*k2)/(-(k1*k2)(1/2)+v0*k2)*k3/(k1*k2)(1/2)/k3)+1)+1/2*k3*log(k1*k2)(1/2)+v0*k2)/(-(k1*k2)(1/2)+v0*k2)*(k1-2*(k1*k2)(1/2)*v0+v02*k2)/k1)/k2

27、9;)g=finverse(f)g=(k1*k2)(1/2)*(exp(-1/2*(2*t*(k1*k2)(1/2)+2*v0*k2+k3*RootOf(exp(-2*(k3*_Z+t*(k1*k2)(1/2)+3*v0*k2)/k3)+2*exp(-(k3*_Z+2*t*(k1*k2)(1/2)+4*v0*k2)/k3)+exp(-2*(t*(k1*k2)(1/2)+v0*k2)/k3)-exp(-1/2*(2*t*(k1*k2)(1/2)+2*v0*k2+k3*_Z)/k3)4-2*exp(-1/2*(2*t*(k1*k2)(1/2)+2*v0*k2+k3*_Z)/k3)2-1)/k3)2

28、-1)/k2/(exp(-1/2*(2*t*(k1*k2)(1/2)+2*v0*k2+k3*RootOf(exp(-2*(k3*_Z+t*(k1*k2)(1/2)+3*v0*k2)/k3)+2*exp(-(k3*_Z+2*t*(k1*k2)(1/2)+4*v0*k2)/k3)+exp(-2*(t*(k1*k2)(1/2)+v0*k2)/k3)-exp(-1/2*(2*t*(k1*k2)(1/2)+2*v0*k2+k3*_Z)/k3)4-2*exp(-1/2*(2*t*(k1*k2)(1/2)+2*v0*k2+k3*_Z)/k3)2-1)/k3)2+1)x,y=dsolve('Dx=0.

29、0439-1.4928*x2,Dy=x','x(0)=v0,y(0)=0')x,y=dsolve('Dx=1-x2,Dy=x','x(0)=v0,y(0)=0')附录二：问题二SAS数据分析程序data jiedianliuliang;input x y;d=x-y;cards;12 6.26 76.53 122.7168.82 75.52 16.3 17.65 14.33 15.51 12 11.39 16.9 16.08 18.06 21.0669 79.0175 73.39 42.08 49.01 17.9 16.8 28.3 23.41 15.3 13.24 18.7 20.64 32.9 32.3127.72 28.89 55.99 56.7825.6 23.48 66.38 63.03 123.2 154.53 14.6 15.37 125.58 152.46 ;run;proc univariate data=jiedianli