一种河道与泛区凌汛壅水-溃堤-淹没耦合模拟方法与流程

一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法
技术领域
1.本发明涉及一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法。


背景技术：

2.极端天气条件下凌汛溃堤灾害是寒区河道冬春季节最突出的重大自然灾害之一，具有险点多、险段长、影响范围广等特点。目前，已有学者开展了溃堤凌汛洪水淹没风险分析二维水动力模型、河道冰塞壅水淹没风险评估一维水动力模型、河冰生消演变数值模拟模型等研究，研究成果只体现在河道、
3.堤防或泛区的单一维度，尚缺少凌汛洪水在河道、堤防、泛区三个维度系统中的耦合模拟方法研究成果，而实际凌汛致灾是河道凌汛壅水、堤防溃决与泛区凌汛洪水淹没的动态耦合过程，目前却鲜有相关研究成果的报道。


技术实现要素：

4.本发明所要解决的技术问题是提供一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法。
5.本发明一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法的技术方案是这样实现的：一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法，依次包括如下步骤：
6.1、数据资料收集整编
7.收集研究区域包括行政区划、居民地、河流水系、交通道路、堤防在内的基础地理数据、历史凌汛洪水过程及灾害资料、数字高程模型即dem、遥感影像等数据资料，并根据凌汛壅水
‑
溃堤
‑
淹没耦合模型构建对数据的要求，对其进行系统化整编处理；
8.2、凌汛壅水
‑
溃堤
‑
淹没耦合模拟方案拟定
9.设定当河道水位达防凌保证水位时堤防溃决，溃口为矩形、瞬间溃决到底，溃口宽度即溃口底宽b为历史实际发生的最大值，分流72h后开始复堵截流，封堵时间24h，选取历史典型年凌汛期溃口上游洪水过程，以研究河段历史最大凌峰流量为峰值进行同倍比放大，将其作为溃堤分析方案的上游入流条件，下游出流条件为凌汛期河道断面的水位
‑
流量关系，溃口下游设置长度10km的全河宽冰塞，凌汛期冰塞糙率n
i
按流量为3200m3/s时凌汛冰塞壅水高度h为2.3m进行设置；
10.3、凌汛壅水
‑
溃堤
‑
淹没耦合模拟模型构建
11.根据河道一维水动力模型和泛区二维水动力模型原理，考虑凌汛冰塞壅水特性及溃堤凌汛洪水演进特征，采用河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法，实现河道凌汛壅水、堤防溃决分流、泛区凌汛洪水淹没的动态耦合数值仿真；
12.(1)、河道一维凌汛壅水模拟模型原理
13.凌汛期河道冰塞造成过水断面河床湿周χ
b
增大、水力半径r减小、河道综合糙率n增加，河道综合糙率n由河床糙率n
b
与凌汛期冰塞糙率n
i
共同决定，由于凌汛期冰塞糙率n
i
与冰塞体特征、凌汛洪水条件、气温等多因素密切相关，并于生消演变过程中动态变化，量
测难度大，故本发明一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法采用凌汛期冰塞糙率n
i
与河道综合糙率n对壅水作用进行优化处理，通过凌汛冰塞壅水高度h对综合糙率n的敏感性分析，确定凌汛期冰塞糙率n
i
以及河道综合糙率n。
14.本发明一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法采用一维非恒定流水动力模型模拟河道凌汛壅水过程，计算不同时刻凌汛冰塞壅水水面线以及不同位置凌汛洪水过流能力，描述河道一维凌汛洪水演进的方程组公式1至公式5依次表示如下：
[0015][0016][0017][0018][0019]
对于天然河道，一般χ
b
≈χ
i
，则：
[0020][0021]
式中：q为河道流量(m3/s)；a为过水面积(m2)；x为沿程距离(m)；t为时间(s)；q为区间源汇项单宽流量(m2/s)；r为水力半径(m)；c为谢才系数(s/m
1/3
)；z为凌汛洪水位(简称水位)(m)，为凌汛冰塞壅水高度(简称壅水高度)h与畅流期(无冰塞)洪水位之和；n为河道综合糙率；n
b
为河床糙率；n
i
为凌汛期冰塞糙率(简称冰塞糙率)；χ
b
为河床湿周(m)；χ
i
为冰塞湿周(m)；α为动量修正系数；v为流速(m/s)；j为水力坡度
[0022]
(2)、凌汛壅水作用下溃堤分流流量计算方法
[0023]
凌汛壅水
‑
溃堤
‑
淹没是河道、堤防与泛区三者之间的凌汛洪水时空动态耦合演进过程，本发明一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法通过溃堤侧向建筑物链接方式实现河道一维凌汛壅水计算模型和洪泛区二维凌汛洪水演进数值模型的实时动态耦合，不同时刻凌汛溃堤分流流量的计算公式(依次为公式6至公式8)如下所示：
[0024][0025][0026][0027]
式中：q
l
为溃口分流流量(m3/s)；b为溃口底宽(m)；h为溃口内侧河道水位(m)；h
b
为溃口底高程(m)；s为溃口边坡比；c
w
为溃口水平部分堰系数；c
s
为溃口边坡部分堰系数；c
v
为入流收缩损失修正系数；k
s
为淹没修正系数；c
b
为无量纲系数；w
r
为溃口处河道宽度(m)；h
d
为溃口最终底高程(m)；q
p
为上一迭代溃口流量(m3/s)；h
ds
为溃口外侧泛区水位(m)；g为重
力加速度(m/s2)；一般c
w
取0.5464，c
s
取0.4319，c
b
取0.7403
[0028]
(3)、泛区二维凌汛溃堤洪水淹没模拟模型原理
[0029]
凌汛溃堤水流湍急，基本不涉及冰盖下水流运动，而且泛区通常具有范围广、地形平坦等特点，因此本发明一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法根据宽浅区域水流挟冰运动特性，利用流凌糙率n
i
与地表糙率n
b
分别计算冰水交界面水流拖拽力与地表摩阻力，数值模拟时采用流凌与地表的综合糙率n对其进行优化处理，凌汛溃堤洪水运动方程公式9至公式14依次表示如下：
[0030][0031][0032][0033][0034][0035][0036]
式中：h为水深值(m)；z
b
为地面高程(m)；u和v分别为x、y方向上的流速分量(m/s)；t为时间(s)；τ
ix
和τ
iy
分别为冰水交界面水流拖拽力在x、y方向的分量；τ
bx
和τ
by
分别为泛区地表摩阻力在x、y方向的分量；n
b
为泛区地表糙率；n
i
为流凌糙率；n为泛区综合糙率，即流凌与地表的综合糙率；ρ为水密度(kg/m3)；g为重力加速度(m/s2)；q
l
为源汇项
[0037]
根据复杂地形地貌条件下凌汛溃堤洪水演进特征和建筑物受灾特点，本发明一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法采用加大糙率法概化房屋建筑物，将阻水效果明显的线状地物作为模型内边界，并对泛区其它土地利用类型进行糙率分区处理，以模拟不同地物与凌汛洪水演进的相互影响。根据干水深(h
dry
)和湿水深(h
wet
)理论，利用淹没水深(h)大小优化网格计算属性，从而提高模型计算效率与稳定性，判断准则为：当h＜h
dry
时，网格不参与计算，质量通量与动量通量均为0；当h
dry
＜h＜h
wet
时，只计算网格质量通量；当h>h
wet
时，同时计算质量通量与动量通量。
[0038]
4、凌汛壅水
‑
溃堤
‑
淹没过程耦合模拟
[0039]
(1)、确定河道凌汛冰塞壅水高度与河道综合糙率的关系
[0040]
选择典型险工段，考虑不同流量条件，分别对应于险工段下游10km河段设定不同凌汛期冰塞糙率n
i
为0.06、0.10、0.12和0.15，基于已构建的河道一维凌汛壅水模拟模型，模拟不同流量条件、不同凌汛期冰塞糙率n
i
组合工况下的凌汛洪水位z及凌汛冰塞壅水高度h变化过程，从而确定河道凌汛冰塞壅水高度h与河道综合糙率n的关系；
[0041]
(2)、凌汛壅水作用下溃堤流量过程计算
[0042]
根据步骤3构建的凌汛壅水
‑
溃堤
‑
淹没耦合模拟模型，设置溃口瞬间溃决到底并于72h后开始封堵、24h后溃口全部合龙，利用嵌入的凌汛壅水作用下溃堤分流流量计算方法，计算溃口宽度及溃口分流流量q
l
变化过程；
[0043]
(3)、泛区凌汛溃堤洪水动态淹没过程模拟
[0044]
基于已构建的凌汛壅水
‑
溃堤
‑
淹没耦合模拟模型，模拟不同方案对应的凌汛壅水
‑
溃堤
‑
淹没动态演进过程，并统计不同时刻凌汛洪水淹没水深h分布情况；
[0045]
(4)、泛区凌汛溃堤洪水淹没风险分析
[0046]
根据凌汛洪水动态淹没过程模拟结果，计算不同时刻凌汛洪水淹没面积以及不同水深等级对应淹没面积变化情况，分析随着凌汛洪水演进时间的推移、淹没面积及淹没范围的变化情况，分析淹没风险较大的区域，明确防凌减灾重点保护范围。
[0047]
本发明一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法具有适于宽浅区域凌汛洪水演进模拟、实现了河道凌汛壅水、堤防溃决分流和泛区凌汛洪水淹没的实时动态耦合、可为凌汛壅水作用下的溃堤淹没风险评估提供重要技术手段的特点。
附图说明
[0048]
下面结合附图和实施例对本发明做进一步的说明。
[0049]
图1是本发明一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法的流程图；
[0050]
图2是研究区域及溃口位置分布图；
[0051]
图3是研究区域地形图；
[0052]
图4
‑
1是凌汛壅水
‑
溃堤
‑
淹没耦合模拟模型入流控制条件示意图；
[0053]
图4
‑
2是凌汛壅水
‑
溃堤
‑
淹没耦合模拟模型出流控制条件示意图；
[0054]
图5
‑
1是q＝2000m3/s情况不同位置凌汛冰塞壅水高度随凌汛期冰塞糙率变化曲线图；
[0055]
图5
‑
2是羊场险工不同凌汛期冰塞糙率情况下凌汛洪水位随流量变化曲线图；
[0056]
图5
‑
3是羊场险工不同凌汛期冰塞糙率情况下凌汛冰塞壅水高度随流量变化曲线图；
[0057]
图6
‑
1是溃口宽度变化过程图；
[0058]
图6
‑
2是溃口分流流量过程曲线图；
[0059]
图7
‑
1是演进时刻12h的凌汛洪水淹没水深分布图；
[0060]
图7
‑
2是演进时刻24h的凌汛洪水淹没水深分布图；
[0061]
图7
‑
3是演进时刻72h的凌汛洪水淹没水深分布图；
[0062]
图8是不同时刻凌汛洪水淹没面积变化曲线图；
[0063]
图9是不同水深等级凌汛洪水淹没面积统计图。
具体实施方式
[0064]
实施例1
[0065]
1、软件来源
[0066]
python：版本3.8；
[0067]
编译器：pycharm
[0068]
2、本发明一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法
[0069]
如图1所示，一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法，依次包括如下步骤：
[0070]
(1)、数据资料收集整编
[0071]
以黄河宁蒙段凌汛灾害发生最为频繁、影响损失最为严重如图2所示的巴彦高勒至头道拐(简称巴
‑
头)河段为研究对象，该河段北岸堤防防洪标准为50年一遇，而南岸堤防防洪标准多为30年一遇，考虑最不利情况，确定以南岸堤防及泛区作为具体研究区域，收集上述具体研究区域包括行政区划、居民地、河流水系、交通道路、堤防在内的基础地理数据、历史凌汛洪水过程及灾害资料、数字高程模型即dem、遥感影像等数据资料，并根据凌汛壅水
‑
溃堤
‑
淹没耦合模型构建对数据的要求，对其进行系统化整编处理；
[0072]
(2)、凌汛壅水
‑
溃堤
‑
淹没耦合模拟方案拟定
[0073]
综合考虑历史溃堤险段、凌汛壅水水位超过防凌保证水位、堤外居民或经济集中等因素，确定堤防溃口位置如图2所示，并基于历史冰塞特征，考虑最不利情况，设定当河道水位达防凌保证水位时堤防溃决，溃口为矩形、瞬间溃决到底，溃口宽度为1990年以来最大的100m，分流72h后开始复堵截流，封堵时间24h，选取2007～2008年度凌汛期奎素段溃口(三湖河口站)上游典型洪水过程，以头道拐站历史最大凌峰流量3500m3/s为峰值进行同倍比放大，将其作为溃堤分析方案的上游入流条件，下游出流条件为头道拐站2007～2008年度凌汛期稳封阶段的水位
‑
流量关系，溃口下游设置长度10km的全河宽冰塞，凌汛期冰塞糙率n
i
按流量为3200m3/s时凌汛冰塞壅水高度h为2.3m进行设置，通过凌汛冰塞壅水高度h对河道综合糙率n的敏感性分析，确定凌汛期冰塞糙率n
i
以及河道综合糙率n。拟定的计算方案列表，见表1
[0074]
表1计算方案列表
[0075][0076]
(3)、凌汛壅水
‑
溃堤
‑
淹没耦合模拟模型构建
[0077]
利用巴
‑
头河段2012年10～11月份实测断面资料，构建河道一维凌汛壅水模拟模型，巴
‑
头河段河长540km，控制断面146个，根据堤防溃口位置确定泛区计算范围，采用非结构化三角形网格剖分区域地形，最大网格边长100m，基于航天飞机雷达地形测绘使命(srtm)90m dem构建不同溃口南岸泛区二维凌汛溃堤洪水淹没模拟模型，并连接建立河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟模型，基础数据统一为中国大地坐标系(cscs)2000高斯
‑
克吕格(gauss
‑
kruger)投影、1985国家高程基准，如图3所示；
[0078]
堤防溃口下游10km全河段凌汛期冰塞糙率n
i
设定为0.15，泛区流凌与地表的综合糙率n为居民地0.08(加大糙率)、旱地0.04，干水深h
dry
为0.005m、湿水深h
wet
为0.1m，一二维
耦合模型计算时间步长为10s，当河道水位达防凌保证水位(堤顶高程
‑
2m)时，堤防瞬间溃决到底，矩形溃口宽度100m，河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟模型的上游入流条件为2007～2008年奎素段溃口(三湖河口站)上游流量过程的同倍比放大数据(凌峰流量为头道拐站历史最大3500m3/s)，下游出流条件为2007～2008年度头道拐站凌汛期稳封阶段的凌汛洪水位
‑
流量关系，如图4
‑
1、4
‑
2所示，相关数据见表2
‑
1、2
‑
2；
[0079]
表2
‑
1 2007～2008年奎素段溃口(三湖河口站)上游流量过程的同倍比放大数据
[0080]
日期时间流量(m3/s)日期时间流量(m3/s)日期时间流量(m3/s)2008/3/1816:001188.822008/3/2008:003308.022008/3/2100:002671.522008/3/1820:001203.592008/3/2009:003393.672008/3/2102:002597.682008/3/1900:001270.042008/3/2010:003359.702008/3/2105:002464.772008/3/1906:001314.352008/3/2011:003305.062008/3/2108:002373.212008/3/1910:001321.732008/3/2012:003325.742008/3/2112:002627.222008/3/1916∶001247.892008/3/2013:003214.982008/3/2200:002265.402008/3/1920:002141.352008/3/2014:003476.372008/3/238:00865.402008/3/1922:002289.032008/3/2015:003337.552008/3/248:00785.652008/3/1923:002348.102008/3/2016:003138.192008/3/258:00843.252008/3/2001:002407.172008/3/2017:002894.512008/3/268:00869.832008/3/2002:303012.662008/3/2018:002761.602008/3/278:00818.142008/3/2004:003046.622008/3/2019:002687.762008/3/288:00776.792008/3/2005:003046.622008/3/2020:002835.442008/3/298:00775.322008/3/2006:003401.052008/3/2021:002791.142008/3/308:00989.452008/3/2007:003325.742008/3/2022:002776.372008/3/318:001083.97
[0081]
表2
‑
2 2007～2008年度头道拐站凌汛期稳封阶段凌汛洪水位
‑
流量关系
[0082]
凌汛洪水位(m)987.45987.76988.00988.25988.53989.29流量(m3/s)280316420450540760
[0083]
(4)、凌汛壅水
‑
溃堤
‑
淹没过程耦合模拟
[0084]
①
、确定河道凌汛冰塞壅水高度与河道综合糙率的关系
[0085]
选择羊场险工作为典型险工段，考虑不同流量q＝400m3/s、800m3/s、1200m3/s、1600m3/s、2000m3/s、2400m3/s、2800m3/s和3200m3/s，分别对应于险工段下游10km河段设定不同凌汛期冰塞糙率n
i
为0.06、0.10、0.12和0.15，组合方案32个，基于巴彦高勒至头道拐段河道一维凌汛壅水模拟模型，模拟不同流量q、不同凌汛期冰塞糙率n
i
组合工况下的凌汛洪水位z及凌汛冰塞壅水高度h变化过程，如图5
‑
1至5
‑
3所示，相关数据见表3
‑
1至3
‑
3，与畅流期相比，随着凌汛期冰塞糙率n
i
增大，冰塞上下游河段均发生不同程度的壅水现象；同一位置凌汛洪水位z随着流量q的增大而升高，相同流量q条件下凌汛期冰塞糙率n
i
越大，凌汛洪水位z与凌汛冰塞壅水高度h越高；相同凌汛期冰塞糙率n
i
条件下，凌汛冰塞壅水高度h随流量q增大而缓慢增加，流量q为3200m3/s，凌汛冰塞壅水高度h为2.3m时对应凌汛期冰塞糙率n
i
约为0.15；由于不同断面的形状各不相同，不同断面凌汛冰塞壅水高度h随凌汛期冰塞糙率n
i
变化的规律也不尽相同；
[0086]
表3
‑
1羊场险工q＝2000m3/s里程与壅水高度曲线数据(n
i
＝0.15)
[0087][0088]
表3
‑
2羊场险工流量与水位曲线数据(n
i
＝0.15)
[0089][0090]
表3
‑
3羊场险工流量与壅水高度曲线数据(n
i
＝0.15)
[0091][0092]
②
、凌汛壅水作用下溃堤流量过程计算
[0093]
根据巴彦高勒至头道拐河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟模型，设置溃口瞬间溃决到底并于72h后开始封堵、24h后溃口全部合龙，溃口宽度变化及溃口分流流量q
l
变化过程，如图6
‑
1与6
‑
2所示，相关数据见表4
‑
1与4
‑2[0094]
表4
‑
1时间与溃口宽度曲线图数据
[0095][0096]
表4
‑
2时间与溃口分流流量曲线图数据
[0097]
(2008年同期)时间溃口分流流量(m3/s)2008/3/2620:0002008/3/2705:00108.502008/3/2710:00106.912008/3/2800:00102.452008/3/2816:0072.482008/3/2819:0067.452008/3/2900:0058.342008/3/2920:0019.452008/3/3006:007.292008/3/3012:000
[0098]
③
、泛区凌汛溃堤洪水动态淹没过程模拟
[0099]
基于已构建的凌汛壅水
‑
溃堤
‑
淹没耦合模拟模型，模拟凌汛壅水
‑
溃堤
‑
淹没动态演进过程，并统计不同时刻凌汛洪水淹没水深分布情况，如图7所示；
[0100]
④
、泛区凌汛溃堤洪水淹没风险分析
[0101]
根据凌汛溃堤洪水动态淹没过程模拟结果，计算不同时刻凌汛洪水淹没面积以及不同水深等级对应淹没面积变化情况，如图8和图9所示,相关数据见表5与表6。分析可知：随着凌汛洪水演进时间的推移，淹没面积不断增大，但增大速率缓慢减小，且基本于75h达
到淹没平衡状态；淹没区不同水深等级对应淹没面积占比分别为(0,0.5]占51.25％、(0.5,1.0]占29.28％、(1.0,2.0]占9.71％、(2.0,3.0]占8.23％、(3.0,5.0]占1.53％，说明泛区地形整体较为平坦、坡度较缓，但凌汛洪水淹没水深大于0.5m的面积占比约为48.75％，可见泛区凌汛溃堤洪水淹没风险较大，淹没水深大于0.5m的近居民区域，应作为防凌减灾重点保护范围。
[0102]
表5洪水演进时刻与淹没面积曲线图数据
[0103]
洪水演进时刻(h)10305080100140淹没面积(km2)4.9313.7318.5220.6420.7720.89
[0104]
表6不同水深等级凌汛洪水淹没面积统计图数据
[0105]
水深等级(m)不同方案对应凌汛洪水淹没面积(km2)＜0.510.710.5～1.06.121.0～2.02.032.0～3.01.72＞3.00.32
[0106]
实施例说明，本发明一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法在收集黄河宁蒙段凌汛灾害发生最为频繁、影响损失最为严重巴彦高勒至头道拐(简称巴
‑
头)河段防洪标准多为30年一遇的南岸堤防历史相关数据并根据凌汛壅水
‑
溃堤
‑
淹没耦合模型构建对数据的要求，对其进行系统化整编处理后，设计了凌汛壅水
‑
溃堤
‑
淹没耦合模拟方案，设定当河道水位达防凌保证水位时堤防溃决，溃口为矩形、瞬间溃决到底，溃口宽度为1990年以来最大的100m，分流72h后开始复堵截流，封堵时间24h，选取2007～2008年度凌汛期奎素段溃口(三湖河口站)上游典型洪水过程，以头道拐站历史最大凌峰流量3500m3/s为峰值进行同倍比放大，将其作为溃堤分析方案的上游入流条件，下游出流条件为头道拐站2007～2008年度凌汛期稳封阶段的水位
‑
流量关系，溃口下游设置长度10km的全河宽冰塞，凌汛期冰塞糙率n
i
按流量q为3200m3/s时凌汛冰塞壅水高度h为2.3m进行设置，根据河道一维水动力模型和泛区二维水动力模型原理，考虑凌汛冰塞壅水特性及溃堤凌汛洪水演进特征，采用河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法，实现河道凌汛壅水、堤防溃决分流、泛区凌汛洪水淹没的动态耦合数值仿真；根据凌汛溃堤洪水动态淹没过程模拟结果，计算不同时刻凌汛洪水淹没面积以及不同水深等级对应淹没面积变化情况，进行泛区凌汛溃堤洪水淹没风险分析后得出泛区凌汛溃堤洪水淹没风险较大，淹没水深大于0.5m的近居民区域，应作为防凌减灾重点保护范围的结论。因而，本发明一种河道与泛区凌汛壅水
‑
溃堤
‑
淹没耦合模拟方法具有适于宽浅区域凌汛洪水演进模拟、实现了河道凌汛壅水、堤防溃决分流和泛区凌汛洪水淹没的实时动态耦合、可为凌汛壅水作用下的溃堤淹没风险评估提供重要技术手段的特点。