研究生数学建模竞赛优秀论文选《唐家山堰塞湖泄洪问题的建模与解决方案》37页.docx
- 文档编号:14645481
- 上传时间:2023-06-25
- 格式:DOCX
- 页数:68
- 大小:547.46KB
研究生数学建模竞赛优秀论文选《唐家山堰塞湖泄洪问题的建模与解决方案》37页.docx
《研究生数学建模竞赛优秀论文选《唐家山堰塞湖泄洪问题的建模与解决方案》37页.docx》由会员分享,可在线阅读,更多相关《研究生数学建模竞赛优秀论文选《唐家山堰塞湖泄洪问题的建模与解决方案》37页.docx(68页珍藏版)》请在冰点文库上搜索。
研究生数学建模竞赛优秀论文选《唐家山堰塞湖泄洪问题的建模与解决方案》37页
全国第五届研究生数学建模竞赛
题目唐家山堰塞湖泄洪问题的建模与解决方案
摘要:
本文通过建立蓄水模型、溃坝模型、演进模型和人员调度模型来分析和解决唐家山堰塞湖的泄洪问题。
在蓄水模型中,假定堰塞湖横截面面积与时间服从分段线性关系,利用微积分的方法得到水位高程和蓄水量之间的函数关系表达式,同时预测出5月25日到6月12日的水位高程,预测值与实测值十分逼近,说明该模型是可行的。
在溃坝模型中,利用微元法推导出逐步溃坝模式下的流量公式,并采用迭代法求出坝址处流量过程线和坝前水位变化趋势。
演进模型在溃坝模型的流量过程线基础上,利用Newton-Raphson迭代法求解圣维南方程组,得到洪水在河道中的演进过程,即洪峰沿河道到达下游各城镇的时间、流量和水位。
根据该演进模型,确定其可能被淹没的区域以及需要撤离的人员数量,建立人员调度的网络流模型。
人员的撤离遵循“步行为主”和“就近安置”的原则,采用动态模拟的方法,得到人员撤退的可行调度方案。
在三分之一瞬时溃坝情况下,运用上述演进模型和人员调度模型,估计出堰塞湖下游共有15个城镇会受到不同程度的淹没,需要撤离的人数约为30万。
根
据上述调度原则,设置8个安置点,建立相应的人员调度网络,并给出相应的撤离方案。
最后对所建立的模型进行了优劣性分析,并给出了模型的改进和扩展方法。
此外,还对地震后各种次生山地灾害的防治以及科研人员的工作提出相应的建议。
关键词:
圣维南方程组;动态模拟;网络流算法;微元法;Newton—Raphson
参赛密码
(由组委会填写)
参赛队号
一、问题的提出
1.1背景
2008年5月12日14:
28在我国四川汶川地区发生了8.0级强烈地震,给人民生命财产和国民经济造成了极大的损失。
地震引发的次生灾害也相当严重,截至5月22日,震灾区共发现堰塞湖33处,相应的位置、坝体特征、上下游城镇分布、灾情状况等详细信息见附录1。
其中以唐家山堰塞湖尤为严重。
图1堰塞湖形成示意图
唐家山堰塞湖距离北川县城6公里,坝顶高程750.2米,坝高82.8米,顺
河长约220米,湖上游集雨面积3550平方公里,地震后每天新增至少500万立方米的水量,同时要承接来自北川和平武方向的来水,随着水位不断上涨,对下游地区造成了严重威胁,见图1。
根据江油市水务局的资料显示,如果唐家山堰塞湖一旦决口,江油市临近北川的6个乡镇和3个重点企业将会受到直接威胁,同时会威胁三个水电站,即通口水电站、香水电站和青莲水电站的正常运行。
为了减轻次生灾害,在有限的时间里挖出泄洪渠道,采取一个简单的泄洪措施,是最可行的措施。
1.2问题
(1)建立唐家山堰塞湖以水位高程为自变量的蓄水量数学模型。
并以该地区天气预报的降雨情况的50%,80%,100%,150%为实际降雨量建立模型预计自5月25日起至6月12日堰塞湖水位每日上升的高度(不计及泄洪);
(2)在合理的假设下,建立堰塞湖蓄水漫顶后在水流作用下发生溃坝的数学模型,包含缺口宽度、深度、水流速度、水量、水位高程,时间等变量。
(3)根据数字地图给出坝体发生溃塌,造成堰塞湖内1/3的蓄水突然下泻时的洪水水流速度及淹没区域,并在此基础上考虑洪水淹没区域中人口密集区域的人员撤离方案。
(4)分析建立的数学模型所采取对策的正确性和改进的可能性。
讨论为应对地震后次生山地灾害(不限堰塞湖),科技工作中应该设法解决的关键问题,并提出有关建议。
二、基本假设
1)从5月12日之后发生的余震假设不会再对唐家山堰塞湖附近的地形、地貌产生太大的改变,在本文的模型中,忽略这些变化引起的影响;
2)堰塞湖水为理想流体,不可压缩,且无粘性;
3)溃坝过程中,在时间T内溃决口的形状形成并且达到稳定;
4)堰塞湖旁侧入流的动量沿干流流向的分量相对较小,可忽略不计;
5)溃口底面高程及溃口宽度在某时段内呈线性变化;
6)溃口具有基本几何形状:
梯形、三角形、矩形;
7)假设灾民在撤退的过程中都愿意服从调度人员的安排;
8)撤离地区的相关信息,包括交通网、居民点、安全区、人口数量等是可得的;
9)淹没区域人员的撤离以群众为主,不考虑留在堰塞湖上参与救援人员的撤离方案。
三、符号说明
t
时间
T
时间区间段的终点
H(t)
时刻t唐家山堰塞湖的水位高程(m)
V(t)
时刻t唐家山堰塞湖的蓄水量(m3)
S(t)
时刻t唐家山堰塞湖的湖面面积(m2)
Di
分别表示堰塞湖的降水量、河流入水、河水损失
vt
时刻t唐家山堰塞湖的平均降水量(mm/s)
ρ
实际降雨量占天气预报降雨情况的百分比
m
可能会遭受洪水淹没危险的城镇数量
n
可供选择的安置点数量
φi
城镇i可能被洪水淹没的区域占该城镇总面积的百分比
Ni
城镇i的总人口数
Cj
安置点j的容量,即可以接受受灾转移群众的数量
Ai
城镇i的总面积(单位:
平方公里)
sij
城镇i到安置点j的直线距离
四、问题分析
本问题是以汶川大地震为背景,考虑唐家山堰塞湖的处理方案以及溃坝后可能造成的损失,并且考虑相应的人员撤离方案。
通过分析可知,本问题可以在建立溃坝模型的基础上,利用洪水在下游河道中的演进状态来确定人员撤离方案及相关应对措施。
在问题一中,要建立以水位高程为自变量的蓄水量数学模型。
根据微积分的知识,建立以时间t为变量的微分方程,并且利用实际统计的数据,通过估计参数的方法,可以得到水位高程和蓄水量之间的关系。
影响蓄水量的因素主要有几个:
降水量、河流入水、河水损失,其中损失包括植物截流、洼坑存蓄、下渗和蒸发等。
鉴于问题要求考虑的时间不足一个月,所以影响堰塞湖蓄水量的主要是降水量和河流的入水。
根据历年的数据记录可以得到河流的入水量,天气预报可以提供降水量的数据,因此可以建立出水位高程和蓄水量的函数关系。
在问题二中,根据唐家山堰塞湖坝体的具体情况知,其漫顶溃坝过程属于土石坝逐步溃坝模式,这时可以利用堰塞湖的水量、水位等相关参数模拟出溃坝过程中的流量、流速等。
问题三在问题二的基础上,首先在三分之一溃坝模式下,由圣维南方程组推算出溃坝洪水在沿途的演进过程,随后根据所得数据确定唐家山下游城镇可能被淹没的程度以及各城镇需要撤离的人员数量,建立人员调度的网络流模型。
通过分析唐家山下游地区地震后道路的情况,以及可能被淹没区域的实际海拔高度,对人员的撤离采取“步行为主”和“就近安置”的原则,人员的财产和物资的撤退以机动车运输为主,以在短时间内完成撤退任务为目标,建立人员撤退的调度方案。
在问题四中,对所建立的蓄水模型、溃坝模型、演进模型和调度模型进行可行性分析并且讨论改进的可能性。
地震后的次生灾害,其危险和危害,甚至可能超过大地震、堰塞湖,因此需要对其它情况进行考虑,比如:
山洪灾害,山体滑坡,泥石流等等。
为了建立完善的数学模型,下面对各个问题分别作详细的分析。
五、模型的建立与求解
5.1问题一的模型建立与算法设计
5.1.1堰塞湖水位高程和蓄水量的关系模型
堰塞湖水位高程,指的是堰塞湖的水面标高,一般以黄海海面为起算点,其水位高程为0米)。
唐家山堰塞湖的横纵剖面图见图2。
图2唐家山堰塞湖横纵剖面图
图2表明唐家山堰塞湖的横纵切面都不是规则的图形,因此在估计水位高程和蓄水量的关系时,不能根据常规的图形进行计算。
横剖面的图形说明唐家山堰塞湖的湖面面积是不断增长的。
设水位高程和蓄水量分别为时间t的函数H(t)和V(t),并且设t时刻唐家山堰塞湖的湖面面积为St。
根据流入堰塞湖中的水体积和堰塞湖新增加的水体积相等,建立出相应的微分方程关系式为:
dV(t)=St⋅dH(t)
(5.1)
根据[4]中的数据可知,唐家山堰塞湖的湖面面积S和时间t之间是分段线性函数,即:
⎧k0t+S(0)0≤t ⎪ S(t)=⎪ (5.2) ⎨ ⎪N-1 (t-TN-1 )+S(TN-1) TN-1 ≤t ⎪⎩kN(t-TN)+S(TN) TN≤t 其中,N表示可测得的时间长度T有限分段的数量,ki(1≤i≤N-1)是常系数。 T0=0,TN=T分别表示可测得的时间起点和终点。 TN到TN+1表示待估计的时间段,取参数kN表示系数ki(1≤i≤N)的平均值,即认为未知的时间段内唐家山堰塞湖的湖面面积S和时间t之间是线性函数。 水位高程H(t)和时间t之间的关系也是一个分段多项式函数,即: H(t)=ai(t-Ti-1 )m+H(T ),Ti-1 ≤t (5.3) i-1 其中H(0)=h0表示在初始时刻堰塞湖的水位高程,ai为常系数。 将St,H(t)的表达式代入公式(5.1),从Ti到Ti+1进行积分,i=1,…,N得: V(Ti+1 )-V(Ti)= Ti+1dV(t)= T Ti+1S(t)⋅dH(t) T ii =maiki⊗Tm+1+aS(T)⊗Tm+1 (5.4) m+1 iiii 其中⊗Ti=Ti+1-Ti为第i+1段的时间间隔。 由前面的分析可知,堰塞湖的水 位上涨主要包含三个因素: 河流基流入水量、降水量和损失,其中损失包括植物 截流、洼坑存蓄、下渗和蒸发等。 设在时间Ti到Ti+1内,三个因素的水量分别为: D1,i,D2,i,D3,i,则堰塞湖水量的增长为: V(Ti+1)-V(Ti)=D1,i+D2,i-D3,i (5.5) 其中,基流入水量DT=q⊗T,q表示河流基流流量,[5]指出 1,iBiB B D2,i-D3,i=0.1ραviS⊗Ti,S表示集雨面积,vi表示第i段的时间间隔内的平均天气预报降雨量,ρ为实际降雨情况占天气预报的降雨量的比例。 其中针对唐家山堰塞湖的具体情况,[4]给出估计的参数q=90m3/s,α值为0.66,S=3550km2,代入表达式(5.5)后,得到: V(Ti+1)-V(Ti)=(qB+0.1ραviS)⊗Ti 将(5.6)代入表达式(5.4),得到: (5.6) mak⊗Tm+1+aS(T)⊗Tm+1=(q +0.1ραvS)⊗T (5.7) m+1 从而得到: iiiiiiBii (m+1)(qB+0.1ραviS) ai= 而 m(k + S(T))⊗Tm (5.8) S(Ti)=ki-1⊗Ti-1+S(Ti-1) (5.9) 给定m后,代入不同的Ti,Ti+1,可以得到相应的ai,i=1,…,N,即得到水位高程H(t)和时间t之间的函数表达式。 对i=1+N,取aN+1为系数ai,i=1,…,N的平均值。 令V(T0)=V0为初始堰塞湖的水容量,对表达式(5.1)从0到t进行积分,设t在第N+1个时间区间内。 N V(t)=V+ aimki(Tm+1-Tm+1)+ tSdH(t') 0∑ i=0 i+1i ⎰Tt' =V+∑ ⎛aiki ⊗Tm+1+S(T )⊗Tm⎞ (5.10) N 0⎜m+1 ii-1i⎟ i=0⎝⎠ +aiki (t-T)m+1+S(T)(t-T)m m+1 iNi 1 ⎛H(t)-h⎞m 根据(5.3),t=⎜0⎟,T ≤t≤T ,代入(7)后得到: ⎝1⎠ N1-1 N1 ⎛⎛⎞m+1⎞ N1-1amk m+1 m+1 aNmkN⎜ H(t)-hm m+1⎟ V(t)=V + ∑ii(T -T)+11⎜⎜0⎟ -T⎟ (5.11) 0 i=0 m+1 i+1i m+1⎜⎝a1⎠⎟ ⎝⎠ 5.1.2蓄水模型的求解 根据参考文献[4]中得到的5月14号到5月22号的回水长度和回水面积数据, 预测从5月25号到6月12号唐家山堰塞湖的水位上涨。 表1唐家山堰塞湖回水长度和回水面积动态变化情况 时间 回水长度(km) 回水面积(km2) 5月14日 7.39 1.14 5月19日 13.51 2.49 5月22日 14.64 3.37 根据式(5.3),坝前水位高程H是关于时间t的分段高阶多项式函数,阶数m不确定。 因此,首先在m取不同的值下估计水位高程,得到一簇对坝前水位高程的估计值,如图3所示,从图中可以看出,随着m的增大,水位高程增加的过程趋缓,且越逼近真实水位上涨的曲线。 当m≥8时,曲线基本不在变化,且最逼近于真实数据,因此,从减小模型的复杂性考虑,阶次m选择为8。 图3m取不同的值情况下坝前水位高程的估计 图4表示根据模型计算出的水位高程和蓄水量随时间的变化图。 可见模型 与真实值相当逼近。 其中时间指标以5月1日为参考点,单位为天,即时间坐标 为40代表6月9日。 值得注意的是,由于真实情况下6月10日开始人工泄洪,因此水位高程与蓄水量急剧下降。 而模型估计时要求不考虑排水渗漏等因素,因此,估计值反映的是不采取任何人工干预的自然情况下堰塞湖水位和蓄水量的增长趋势,该模型可作为堰塞湖自然发展状况的风险预测参考,具有重要的指导意义。 750 x104 3 7402.5 7302 7201.5 7101 700 1020304050 t 0.5 1020304050 t 图4唐家山堰塞湖水位高程和蓄水量随时间变化的关系图 由模型5.1.1估计出的以水位高程为自变量的蓄水量变化图如图5所示。 由于模型不考虑渗漏泄洪等,因此作为参考,真实值的范围仅考虑到水位上升阶段。 由图可见,蓄水量与水位高程的关系与真实的情况非常接近,说明该模型具有很高的估计精度。 x104 3 2.5 2 1.5 1 0.5 700705710715720725730735740745750 水位高程 图5ρ=100%时唐家山堰塞湖水位高程估计效果图 根据资料,5月14日到6月12日的天气预报降雨量如附录3所示。 当真实降雨量分别为预报降雨量的50%,100%,150%,200%时,坝前水位上升高度的估计值如附录4所示,图6反映了坝前水位估计值随着不同的ρ进行变化的情况。 760 740 720 700 ρ=50% 760 740 720 700 ρ=100% 680 1020304050 t ρ=150% 760 740 720 700 680 1020304050 t 680 1020304050 t ρ=200% 760 740 720 700 680 1020304050 t 图6ρ=50%,100%,150%,200%时唐家山堰塞湖水位高程的计算估计值 5.2问题二的模型建立与算法设计 5.2.1逐步溃坝模型 这里采用美国国家气象局编制的溃坝洪水预报DAMBRK模型。 该模型由三部分组成: 1)大坝溃口形态描述。 用于确定大坝溃口形态随时间的变化,包括溃口底宽、溃口顶宽、溃口边坡及溃决历时; 2)水库下泄流量的计算; 3)溃口下泄流量向下游的演进; 溃口是大坝失事时形成的缺口,溃口的形态主要与坝型和筑坝材料有关。 目前对于实际溃坝机理仍不是很清楚,因此溃口形态主要通过近似假定来确定。 考虑到模型的直观性、通用性和适应性,一般假定溃口底宽从一点开始,在溃决历时内,按线性比率扩大,直至形成最终底宽。 若溃决历时小于10分钟,则溃口底部不是从一点开始,而是由冲蚀直接形成最终底宽。 溃口形态描述主要由四个参数确定: 溃决历时、溃口最终底宽、溃口边坡、溃口底部高程。 由第一个参数可以确定大坝溃决是瞬溃还是渐溃。 由后面三个参数可以确定溃口断面形态为矩形、三角形或梯形及局部溃或全溃。 1.唐家山堰塞湖基本情况 根据现场查勘资料,唐家山堰塞湖基本情况如下: 唐家山堰塞湖坝体位于北川县城上游6km的通口河上,下距苦竹坝约1km,集雨面积3550km2,堰塞湖总容积约3.15亿m3。 坝址区通口河为不对称的V型谷,右岸较陡,坡度60º左右;左岸坡较缓,坡度约为30°。 坝区两岸基岩为寒武系下统清平组薄层硅质岩、砂岩、泥灰岩、泥岩组成。 岩层软硬相间,产状N60°E/NW∠60°。 左岸为逆向坡,右岸为顺向坡。 基岩裂隙较发育,岩体较破碎,强风化带厚度5~10m,其下为弱风化,高程大体为720m。 两岸坡分布有残坡积的碎石土,厚度0~15m。 碎石土由粉质壤土、岩屑和块石组成,其中粉质壤土占60%左右,岩屑占30%~35%,块石占5%~10%。 2.溃口流量公式的推导 土石坝、拱坝、混凝土重力坝的溃决形式以及溃决过程是不一样的。 拱坝的溃决过程较快,几乎在瞬间发生,故模型按瞬溃考虑;混凝土重力坝则可能有一个或者多个坝段溃决,因此其溃决时间稍长,约为几分钟;土石坝不会全部溃决,也不会瞬间溃决,溃决时间相对最长,约从几分钟到几个小时,溃决的快慢,取决于土石坝的土力学参数。 这里由唐家山堰塞湖坝体的成分可知,堰塞湖坝体为土石坝。 在溃决过程中和溃决后的溃口的形状不同,其流量和水位的过程也不同,因此有必要分成两个过程分别进行研究。 图7为溃口大体形状。 图7大坝发生溃决的溃口截面图 溃决过程中: 在发生地震时间T内溃决口的形状形成并且达到稳定,我们假设坝的高度和溃口宽度随着时间呈线性变化,即Hd=k2t,Wup=k1t又因为Wdw=kWup,所以 Wup+Wdw (1+k)⋅k⋅k ⋅t2 A=Hd=12 22 ,又v=Cv ,由此便得到了在发生地震 溃决口的形状形成的时间t内的流量与水位的关系: Q=A⨯v=(1+k)kkt2⨯C (5.12) 全部溃决后: 212v 根据流体力学中的基本知识,在全部溃决后,其缺口面积基本保持不变,由此建立贝努力流体方程 pv2pv2 γ0+0+H=1+1 (5.13) 2gγ2g 其中 p0: 水库的液面上的大气压; p1: 大坝缺口处的侧面液面大气压强; v0: 水库的水流速; v1: 大坝缺口处的水流速; H: 水库与大坝缺口处的高度差; 因为水库的水可以看作静止的,所以v0=0,水库的大气压与大坝缺口处的大气压相等,p0=p1。 将上述条件代入式(5.12)可以得到洪水流出速度为: v1=。 由方程(5.13)推出的洪水流出速度v1,是在理想情况下得到的。 在现实生活中,由于各种因素影响,我们必须在公式中加入系数修整,令Cv为行进流速的 修正,通常为0.96∼0.99之间,即v1=Cv。 根据溃口处洪水流量与洪水流出速度及溃口宽度之间的物理关系,建立微分 方程模型: ⊗Qb=⊗W⊗v 于是得到时刻t的溃决口流量为: (5.14) Q(t)= H(t) ⎰ Hb W(H)v(H)dH= H(t) ⎰ Hb W(H)Cv 2gHdH (5.15) 模型的简化处理: 由于土石坝发生溃决其溃口有不同的形状,可能是三角形、梯形,也可能是矩形(我们对溃口进行了简化,限定其为几种基本的几何图形),以下通过不同的缺口模型计算其流量。 (a) (b) 图8土石坝发生溃决后溃口的形状 (c) 为了简化计算,我们用参数K来表示溃口的上底宽与堤口底宽的比值,如图 8,根据上述流量的物理关系模型及图形的几何关系,可以设计以下理想模型: W(H)=Wdown+(Wup-Wdw)⨯(H-Hd)Hd (5.16) 令Wup=k⨯Wdown,则W(H)=Wdown+(k-1)⨯Wdown⨯(H-Hd)Hd,可得流量 H(t)⎧ 11⎫ Q(t)= ⎰⎨WdownCv 2gH2+(k-1)WdownCv2gH2(H-Hb)/Hb⎬dH (5.17) Hb⎩⎭ 当k→∞时,缺口模型为一个三角形,如图8(a)所示;当1 示。 引入具体的流体力学中的参数: Ks: 尾水影响出流的淹没修正;bi: 溃口的瞬时底宽;z: 溃口边坡比。 则积分(5.17)变为 1.52.5 Q(t)=KsCv(3.1bi(H-Hb) 这就是溃口流量模型。 5.2.2逐步溃坝模型的求解 +2.45z(H-Hb)) (5.18) 对于逐步溃坝而言,需要确定两个参数: 最终溃口底宽b和形状参数z。 最终溃口底宽b和平均宽度b关系为b=b-zhd。 假设溃口底端从一点开始,在溃决时间T内以线性速度增长,直到形成最终溃口b和最终溃口底高程hbm。 这里hb和bi满足: h=h -(h- tρ=tρ (5.19) bdd hbm)(T),bi b() T 其中0≤t≤T,hbm表示最终溃口底高程,ρ为非线性程度参数。 当ht-hb h-hb >0.67时 k=1.0-27.8⎛ht-hb-0.67⎞ (5.20) s⎜h-h⎟ ⎝b⎠ 否则ks=0,式中ht为尾水水位高度。 Cv由下式计算: Q2 Cv=1.0+0.023b (5.21) B2(h-h)2(h-h) dbmb 这里Bb为水库坝前宽度。 水库水量动态平衡方程: ⊗V
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 唐家山堰塞湖泄洪问题的建模与解决方案 研究生 数学 建模 竞赛 优秀论文 唐家山堰塞湖 泄洪 问题 解决方案 37