数学建模 美赛特等奖论文(中文版)分析溃坝:针对南卡罗来纳州大坝坍塌建立模型 下载本文

分析溃坝:针对南卡罗来纳州大坝坍塌建立模型

摘要

萨鲁达大坝建立在卡罗莱纳州的墨累湖与萨鲁达河之间,如果发生地震大坝就会坍塌。

本文通过建立模型来分析以下四种大坝决口时水的流量以及洪水泛滥时水的流量:

? 大坝的绝大部分被瞬间侵蚀看成是大坝瞬间彻底坍塌; ? 大坝的绝大部分被缓慢侵蚀看成是大坝延期彻底坍塌; ? 管涌就是先形成一个小孔,最终形成一个裂口; ? 溢出就是大坝被侵蚀后,形成一个梯形的裂口。

本文建立了两个模型来描述下游洪水的泛滥情况。两个模型都采用离散网格的方法,将一个地区看成是一个网格,每个网格都包含洪水的深度和体积。复力模型运用了网格的速度、重力以及邻近网格的压力来模拟水流。下坡模型假定水流速度与邻近网格间水位高度的成正比例。

下坡模型是高效率的、直观的、灵活的,可以适用于已知海拔的任何地区。它的两个参数稳定并限制了水流,但该模型的预测很少依赖于它们的静态值。

2对于萨鲁达溃坝,洪水总面积为106.5km;它还没有到达国会大厦。罗威克里克的洪水向上游延伸了4.4km,覆盖面积达1.6?2.4km2

变量及假设

表1说明了用来描述和模拟模型的变量,表2列出了模拟程序中的参数。

表 1模型中的变量.

变量 定义 溃坝时的水流量速率

QTF1 瞬间彻底坍塌 QTF2 延期彻底坍塌 QPIPE 管涌 QOT 溢出 Qpeak 最大流速 溃坝时水流出到停止所用时间

tTF1 瞬间彻底坍塌 tTF2 延期彻底坍塌 tPIPE 管涌 tOT 溢出

?V 溃坝后从墨累湖里流出的水的总体积 VolLm 墨累湖的原来体积 AreaLM 墨累湖的原来面积 dbreach 从裂口到坝顶距离

tbreach 从裂口开始到溃坝形成的时间 m 近似圆锥的墨累湖的侧面

1

一般假设

? 正常水位是在溃坝前的湖水位置。 ? 河道中的水流不随季节变化而变动。

? 墨累湖里的水的容积可以看作为一个正圆锥(图1 )。

表2 模拟程序中的参数

参数 所取值 意义

BREACH_TYPE 变量 瞬间彻底坍塌,延期彻底坍,管涌,溢出模型中的一种

?T 10.0 时间不长的长度(s)

MIN_DEPTH 0.0001 网格空时的水的深度(m) TFINAT 100000 大坝彻底决口所用时间 Tb 3600 溃坝达最大值的时间(s) Qpeak 25000 溃坝的最大流速(m3/s) dbreach 30 蓄水池的最初深度(m)

9VolumeLM 2.714?10 墨累湖的总体积(m3)

AreaLM 202?106 墨累湖的总面积(m2)

k 0.504 扩散因素 (控制两网格间交换的水的数量) MAX_LOSS_FRAC 0.25 单位网格中水的最大流失量

图 1. 水库近似一个正圆锥.

大坝假设

? 萨鲁达大坝在以下四种方式之一坍塌: -瞬间彻底坍塌, -延期彻底坍塌, -管涌, -溢出。

2

? 土建大坝的材料是统一的。

? 裂口的基本宽度是在大坝高度和3倍大坝高度之间[美国美国陆军工程1997]

? 没有什方法能防止溃坝。

下游假设

? ? ? ? ? ?

忽略建筑物对水流的阻力。 洪水对地形的影响不显著。

洪水流经的地方不产生任何冲击物。 忽略洪灾前该流域的水。 只考虑溃坝时流出的水。

没有什么方法能防止洪水灾害。

已知数据

? ? ? ? ?

墨累湖的面积:200km2

墨累湖的体积:2.710?106m3

大坝的高:63.4m (坝顶高于海平面370ft) 大坝的长: 2.4km

墨累湖的海拔:高出海平面106.5?110m

模型设计

溃坝

各类溃坝是以流量来描述的,它是时间和对应参数的函数。 瞬间彻底坍塌

瞬间彻底坍塌的流量模型是正三角形[美国美国陆军工程1997](图2)。参数是裂口深度和流出最高量,其值为

dbreach?20m,Qpeak?30,000m3s.

延迟彻底坍塌

利用等腰三角形模型来描述延迟坍塌是可行的,因为它把水量总体积的一半从被侵蚀的大坝中流出直至决口水的流速才达到最大值 [美国美国陆军工程1997] (图3 )。此外,土建大坝的侵蚀时间可能会比其他类型的大坝更长,如混凝土大坝。 这种模型有同样的参数,并有相同的值:

dbreach?20m,Qpeak?30,000m3s.

3

图2瞬间彻底坍塌的流速

图3延迟彻底坍塌的流速

管道崩裂

对于管道崩裂,裂口首先出现在坝面中部,继续增大直至管道上的材料完全瓦解[泥沙淤积和河道水力学组2004年]。由于裂口增大,流速呈指数增加;当管道上的材料完全瓦解时出现最高水流速度。从这个角度看,溃坝时水从流出大坝到静止类似一个大坝完全崩裂的过程。我们利用一个递减的指数函数从彻底崩裂模型中得到不同的结果(图4)。

我们利用增长速度来确定洪峰发生的裂口时间,由于洪水慢慢退去最终水流速度小于洪峰流速的1%。将裂口深度、大坝的洪峰流出量以及裂口时间作为参数,其值为:

dbreach?20m,Qpeak?30,000m3s,tbreach?50,000s.

4

为了更好说明当裂口开始形成时流速与时间的关系,我们绘制出短期内的速度变化情况,如下图4:。

过难关

图4管涌崩塌时的流速

图5管涌崩塌开始时的流速

溢出崩塌

对于溢出崩塌,水开始从裂口的顶部流过,就是说从上面侵蚀着大坝。我们找到关于溢出崩塌的资料不多。在管涌失败中,根据抛物线的形状,我们估计流速增加,直至大坝完全被侵蚀(图6)。在到达裂口时间后,就认为流量等于完全崩塌状态时的大小。

参数仍然是裂口深度,大坝流出量的峰值,以及裂口时间,其值为:

dbreach?20m,Qpeak?30,000m3s,tbreach?30,000s.

5

图6溢出失败的流速

下游流量

我们根据下游地区的水的流动规律建立模型,使用离散的方法建立裂口模型。根据伯努利方程,复力模型采用了物理类比方法检测水的流动;下坡模型对水流量采用了更简单、更直观的水流分析机理。因为复力模型的结果不符合常理;因此,我们在分析洪水灾害时利用下坡模型。

对于这两个模型,将围绕萨鲁达大坝的区域划分成一系列正方形的网格单元。每个单元涵盖面积为210m?210m,并和海拔高度和水的体积(基于单元中水的深度)有关系。该海拔数据摘自美国地质调查所的国家海拔数据[2004](以减少过程时间),平均每组7?7个元胞。每个模型以单元模拟了之间水的泛滥;模型区别出每段时间内相邻元胞之间水是如何交换的以及交换量是多少。

复力模型

构建模型

这个模型对单元内所含的水进行力学分析。每个单元有一个相关的海拔、水的平均深度、平均流速(x?y组成)。外力作用于某一特定单元,假定只产生两种影响:外界压力的施加使四个元胞直接接触,再加上引力的作用,加快了水流向海拔较低的地方(即是下坡) 。 模型的主要原则是:

? 单元之间的流量是与四个相邻共面单元之间不同的压力成比例的; ? 单元间的压强是与平均深度成比例的。

如图7,8所示,把单元压缩至一半深度所产生的压强假设为平均压强,记

1做P??gd

2

6

图7施加于网格单元上的压力示意图

图8施加于网格单元上的重力示意图

我们假定所施加的压强是两个单元压缩的深度平均后的压强,由于边界的压缩深度取平均值,所以使它们之间水深度的变化是呈线性的,施加在相邻单元的压力就是平均压强乘以两个元胞之间的面积。我们用单元内水的体积乘以密度求出质量,再用水所受到的力除以它的质量,即可求出水的加速度:

22g(d0x?d2x)ax?.4?d1

我们估计出当前及其4个相邻的单元倾斜度进而计算出加速度。从几何上确定水平加速度(图9)是

7

2?h?g.4?2??h2

在模型中包括大量的时间段,通常以1秒钟为单位。在每段时间的开始,水流入到包括溃坝的单元中;水量取决于上述的溃坝模型。对于每个时间段,加速度(x?y组成)是用来计算地区内每个单元内增加的水速,则水速是以下面的式子进行变化的,即

vnew?vold?a?t.

ag?

图9 确定沿斜面向下的重力值

水流的速度将决定每一个单元的流动方向:当vx?0时,单元内水向右流动; 当vx?0时,元胞向左流动。在那个方向上流出的水量和速度是成比例的,所以,当前网格单元水的深度变化公式为:

davg?tddonated?.2?

8

流出水量使相邻单元水的深度也发生了变化,所以模型中水的总量是守恒的(忽略流向边缘和注入崩裂的单元的水)。

当速度非常大时,一个单元可以流出的水比它原来拥有的还多。具体来说,如果速度乘以时间大于单元宽度,则流出量将大于网格单元当前的容积。如果出现这种情况,假设网格元胞失去全部的水分,根据x和y方向上标注的测出失水量可以说明这种情况。 评价

模型特点:它用一个简捷而富有意义的物理类推法为水的状态建立了模型。用到的力学分析包括采用伯努利方程梯度和模拟离散流体。

用计算机模拟并保存速度信息可以模拟出受灾地区水的态势。例如,当水袭击哥伦比亚的一个特别建筑物时,如国务院,该模型可以预测出水的速度。 模型的缺点

结果显示该模型是不现实的。因为含有大量水的单位格加速度小,这些单位格流尽是非常缓慢的,即使相邻的是完全的空单元;同理,小单位格内水流尽太快。结果是一个棋盘格局:大单位格变得更大而它们相邻的小单元格则变得更小。这个错误与我们的假设有关:一个单位格内所有的水有相同的速度:一个单独的单元格内的水不能向各个方向流。模型适用于比较简单的地形(如一个简单的下坡通道);但是对需要在几个方向都传输水这种非常复杂地形就不适用了。

这个模型的另一个缺点是其太复杂。对于每个单位格,该模型都篡改了大量的参数,使调试和故障排除变困难了。 下坡模型 构建模型

下坡模型假定两单元间的流速是与这些单元格的质心高度差和接触面积的乘积成正比。该模型允许单元格的水可多方向地流动,如果它高于相邻单元格。在力学模型中,我们流入溃坝地区的水是在每隔一段时间内增加的。

每个时间段,每一个单位格(除了那些在地图底部和右界线上,这是处理后的)与这两个单位格交换水后将立刻下降。这样可以确保每一个单位格与它的四个邻居精确地交换水。为了与单位格交换水,单位格应按下面的公式改变其高度。

ddonated?kdavg?h0?h2? 假设常数k,在坍塌时水的的速度为30m/s。然后我们运用k的变化求出了模型的结果。

然后,把邻近单位的高度变成相反数。为了确保结果一致性,直到时间段的最后变化了的高度才会被使用,毕竟已经被计算出了变化了的单元的海拔。如果一个单位格的失水量已经超过它的最大失水量,我们就取这个值为最大失水量。最大失水量是用来防止水来会无规律流动:即当大单位内水流向了相邻的单元格而自身含水量就越来越少,到达一定程度后水就自然往回流,这就使在任何时刻单元格的一半总是空的。

对于处在边界的单元,根据假设边界单元的失水量等于地图上相对同样位置上单元的失水量。因为这些单元是远离裂口和重要的区域,所以它们精确度不太重要的。我们的方法确保水到达区域边缘时是平稳的而不是非正常的澎湃状态。

9

评价

该模型计算简捷,通俗易懂,容易检测。尽管方程不能用直接物理类推法来推出单位格的水交换,但它的结果还是满足期望的。水快速地流到山下,或穿过附近的平坦的平原灾区,只有水位上升时才能到达山顶。

检验和结果

检验

为了测试我们的模型,我们从U.S.GeologicalSurvey[USGS2004]引用国际海拔数据。这数据是一系列海拔数据,把它划分成行和列。每个单元表示一个长度为30m的正方形。为了减少计算,我们把单元平均按7?7分组,使我们用的这些单元是边长为210m的正方形。 我们设裂口处在坝面的正前部位,模拟在不同的溃坝方式(瞬间彻底坍塌, 延

,6240s)内洪水的蔓延方期彻底坍塌, 管涌, 溢出)下不同时间段(360,1800式,来同时检验静力模型和下坡模型。我们利用瞬间彻底坍塌模型检测k值的影响因素和最大失水量的影响因素。

为了避免在考虑体积容量较小的单元格时得到不精确的结果,当它的平均深度小于0.0001m我们就把它看成空单元格。这个做法在静力模型中特别重要,因为当这些小单元格与含有大量的水的大单元格相邻时,水流向这些小单元格时会需要极大的水流速度(?106m/s)。

为了充分检验下坡模型,我们将瞬时崩塌模型在50000个时段内分别运行。该模型在超过40000个时段后计算得到流出区域时的水速比大坝崩塌时的水的速度还大(速度应该是减小了),所以在超过40000个时段后模型效果很差。我

们也在大速度的情况下测试了该模型。当单位格内裂口的高度在每时段内以大于10m的速度(比模型模拟出的任何速度都大于30倍)增长时,模拟结果在持续1000个时间段后就变的不稳定了。 结果 洪灾程度

洪灾程度是主要是依赖于大坝裂口的类型(图10 );不同的裂口类型使洪水的蔓延速度也不同。对于瞬时崩塌模型,洪水泛滥最大范围为106.5km2. 洪灾最严重的地区是十分平坦且宽阔的萨鲁达和坎格瑞山谷。从这些山谷来的水给哥伦比亚造成洪灾的可能性是很小的。所以我们就不考虑从遥远的坎格瑞山谷下来的洪水的影响 ,但我们希望这些水在可以比得上模拟区域内的洪水。 罗尔斯河

在这个地区,发生在罗尔斯河水灾是十分频繁的,但不是在上游。虽然很难确定萨鲁达河泛滥的具体位置及罗尔斯河泛滥开端,我们估计罗尔斯河周围地区的泛滥值在1.6到2.4 km2之间。最远点是在上游4.4km处。 国会大厦

即使是在瞬间彻底坍塌最极端的情况下,洪水也不能达到国会大厦。

10

误差分析,灵敏度和鲁棒性

用因素k来衡量该模型水流出的量。我们假设水的速度在坍塌期间的高峰流量为30m/s。且k值并不影响模型计算;当k是在100内变化时,在1800个时间短后的总淹没面积仅占17 %。

最大失水量是用来防止每个单元流水太多。发洪水的程度和它的值关系不大,我们取0.25描述合理的水分布状态。

优点和缺点

优点

该模型可以针对一个地区进行模拟:只要给出这个地区的海拔数据,再以从溃坝中流出水的速度为变量建立方程,就可以从中看出洪水的情况。 下坡模型是很直观的。它可以根据简单的网格单元之间的交换规则容易地检验出故障并调整。调整就是要能够解决更多的问题和更为极端的水灾情况,这就要求能根据模型解出另外的结果或得出其他无法预料的结果。 该方法是很有效率的;在单位时间段内该地区网格单元数量是线性的。这种效率使它可以模拟出裂口类型多种可能的变化,计算出在短时间段内的水流速速。

模型产生了三套数据网格:0/1静态值能描述那些发生了洪灾的网格单元; 每个网格单元内水的深度可确定各个地区的水灾严重性;而水的深度要加上每个网格单元的海拔。从绘制出的数据集图形就能很显然的看出洪水泛滥程度和严重性。

缺点

此模型的主要缺点是没有对被洪水淹没的地势较低地区状态预测出其发展趋势。或许可以通过引入一个由最大失水量来确定的深度值来弥补这种缺陷。 另一个缺点在于时间刻度,但这个缺点可以通过更细致的分析来改良它。由于k-依赖(唯一一处时间间隔用的很清晰的地方)是一个缺点,模型的时间间隔是不易改变的。时间段可校准通过分析系统的运行模型,如以淡化渠道的泛滥,并确定水的速度。由于时间尺度并不需要分析水灾的程度,我们没有进行校准。

References

Chauhan, Sanjay S., et al. 2004. Do current breach parameter estimation techniques provide reasonable estimates for use in breach modeling? Utah State University and RAC engines& Economists. www.engineering.usu. edu/uwrl/www/faculty/DSB/breachparameters.pdf.Accessed3February 2005.

Chereisnoff , Nicholas P. 1981. Fluid Flow: pumps, pipes , and Channels. England: Butterworth .

Federal Energy Regulatory Commission (FERC). 2002. Saluda Dam remediation: Updated frequently-asked questions and answers. www.ferc.gov/industries/hydropower/safety/asludaqa.pdf. Accessed February 2005

11

FreadD.L.1998. Dam-breach Modeling and flood Routing: A perspective on Present Capabilities and Future Directions. Silver Spring, MD: National Weather Service, Office of Hydrology.

Mayer,L.1987。Catastrophic Flooding. Boston ,MA :Allen & Un win .

Munson ,Bruce R., et al. 2002. Fundamentals of Fluid Mechanics of Fluid Mechanics. 4th ed. New York: John Wiley& Sons.

Sedimentation and River Hydraulics Group. 2004. Comparison between the methods used in MIKE11,FLDWAV1.0,AND HEC-RAS 3.11 to compute flows through a dam breach. U.S .Department of the Interior.

Smith, Alan A. Hydraulic theory: Kinematic flood routing. Alan A. SimthInc. http://www.alanasmith.com/theory-Kinematic-Flood-routing.htm. accessed 3 February 2005.

U.S. Army Corps of Engineering and design-Hydrologic engineering requirements for reservoirs. Department of the Army Publication EM1110-2-1420,Ch.16.

U.S. Geological Survey (USGS). 2004. Seamless Data Distribution System, National Center for Earth Resources Observation and Science. seamless. usgs. gov. Accessed 4 February 2005

Wahl, Tony l.1997. Predicting embankment dam breach parameters—A needs assessment. Denver, CO: U. S. bureau of Reclamation. http://www.usbr. Gov/pmts/hydraulics_lab/twahl/publications.html. Accessed 3 February 2005.

Williams, Garnett P.1978. Hydraulic geometry of river cross sections— Theory of minimum variance. Geological Survey Professional Paper. Washington, DC: U. S. Government Printing Office.

Appendix: Dam Breach Model Equations

For an instant total failure:

?Qpeak(t?tTF1),??QTF1(t)??tTF1?0?Where

2?VtTF2?Qpeakt?tTF1;tTF1?t,

For a delayed total failure:

?2Qpeakt1,t?tTF2;?2?ttf2??2(tTf2?t)Qpeak1QTF2(T)??,tTF2?t?0andt?tTF2?0;t2TF2??0tTF2?t,???

Where

2?VtTF2?Qpeak.

12

.

For a piping breach that turns into a total failure:

?(Qpeak?1)t/tb?1,t?tb;?5(t?tb)?QPIPE(t)??Qpeakexp[],tb?t?0andt?t?0;

t?tb??0t?t,?Where

5(QPIPE(t)??V?tbreach[2?QpeakIn(Qpeak?1)?5?1)]

Qpeak(1?e)For an overtopping breach:

2?Qpeak?15(t2?2ttbreach?tbreach)?10?6,t?tbreach;??Qpeak(t?tOT)QOT??,tbreacht?0andt?tOT?0;t?t?breachOT?0,otherwise,?

Where

32(?V?0.000005tbreach?tbreachQpeak)tOT??tbreach.

Qpeak

13