一种滑坡涌浪计算方法与流程

本发明涉及一种滑坡涌浪计算方法，属于水利工程灾害预测领域。

背景技术：

水利水电建设的长足发展为社会发展做出了重要贡献，但也带来一系列地质环境问题。相比于一般的边坡失稳问题，库岸边坡失稳会对下游的人民生命财产安全构成更加巨大的威胁。滑坡体滑入水库后产生的涌浪有可能翻越坝顶，淹没下游地区，甚至由于涌浪的冲击作用，大坝存在溃决的风险。在国内水利水电工程建设中，不乏水库滑坡失稳破坏并产生涌浪灾害的先例。如1961年3月柘溪水库塘岩光滑坡，165万m³的滑坡体冲入库区后产生了高达21m的涌浪，造成重大损失。随着一系列高坝大库在我国西南地区建立，水库库岸边坡稳定问题越来越吸引人们的注意。在此背景下，数值模拟作为灾害预测的手段之一，就显的非常必要和重要了。

目前的数值模拟基本是首先用能量守恒原理估算滑坡体入水前的滑速，而后依据经验公式计算初始涌浪，最后用计算流体力学的方法求解涌浪的传播过程，实际上只是求解了已知初始条件和边界条件的计算流体力学问题。滑坡体入水产生涌浪的过程是一个复杂的水动力学问题，过程中涉及到三相介质的相互耦合、浸入、波浪破碎等现象，介质间的界面难于捕捉描述。现有的数值方法，包括基于浅水方程的解法，不能反映三相介质的强烈相互作用。主要的界面描述方法有VOF法、Level Set法等，涌浪自由面的描述在计算流体力学中本身就是很大的难题，数值计算的准确性可以得到进一步改善。相对于涌浪产生来说，涌浪的传播过程研究成果已较为丰富。但此类计算的前提是已知首浪的各种要素，如浪高、速度等。从某种意义上说，现有方法是用一个粗估的初始涌浪数据作为初始计算条件，即使涌浪传播过程的模拟再精细，也很难得到合理可信的计算结果。

尽管国内外学者对此已做了大量的研究工作，取得了丰富的研究成果，但基本上是着眼于某个方面的研究，还没有针对从库岸边坡失稳运动到涌浪产生与传播全过程系统的数值模拟。而实际上此类灾害的发生基本是在几十秒到几分钟极短的时间内，库岸边坡从失稳、滑动到涌浪产生再到涌浪传播应该是一个完整连续的过程，每个阶段的准确描述都将对最后灾害预报的可信度产生较大影响。

技术实现要素：

本发明所要解决的技术问题是针对背景技术的缺点，将库岸边坡滑坡体入水到涌浪传播各阶段作为完整统一的过程进行分析研究。采用原始的不可压缩粘性流体Navier-Stokes方程描述水体、空气、滑坡体的相互作用。采用改进的守恒式Level Set方法，准确捕捉流体自由表面，真实模拟涌浪产生过程，为涌浪传播及灾害预测提供合理的前提条件。在此基础上，采用沿水深积分的流体控制方程模拟涌浪传播过程，对涌浪运动轨迹和泛滥区域进行数值预报，评价涌浪灾害。本发明是在已经利用动力时程强度折减法进行土质库岸边坡动力稳定分析，确定了库岸可能的失稳滑动范围的前提下进行的。

本发明为实现上述发明目的采用如下技术方案：

一种滑坡涌浪计算方法，包括如下步骤：

步骤1：根据已有地形、水位和滑坡体形状位置资料，建立三维有限元模型；

步骤2：结合Navier-Stokes方程和改进的守恒式Level Set方法，计算流体各运动参数，捕捉自由面；

步骤3：采用沿水深积分的流体控制方程模拟涌浪传播过程，对涌浪运动轨迹和泛滥区域进行数值预报。

所述步骤2中对土质或松散介质坡体，基于“流体化”的特点，以流体力学基本方程Navier-Stokes方程配以非牛顿流体或拟流体本构关系描述其物理特性和运动过程，由于将滑坡体视作似流体，滑坡体、水体和空气三相介质可用同一本构模型表示：

式中，τ_y0表示流体屈服应力，φ为内摩擦角，I_2D是应变张量D的第二不变量，μ为黏性系数，η为表示流体特性的参量。

在利用Navier-Stokes方程计算出的流场信息的基础上，采用改进的守恒式双重Level Set方法捕捉滑坡体、水体和空气三者的交界面，反映三者间强烈的流固耦合作用，真实模拟滑坡体入水产生涌浪的过程，改进的守恒式双重Level Set方法包含两个部分：指示函数的输运，和指示函数的重新初始化。

在三相流中，各点的物理参数用指示函数H表示，各点的物理参数包括密度或黏性；引入指示函数H₁、H₂，滑坡体、水体和空气三相介质各自区域分别记作Ω₁、Ω₂、Ω₃，滑坡体与水体、水体与空气、滑坡体与空气界面依次用Γ₁₂、Γ₂₃、Γ₁₃表示，P表示物理参数，则各区域参数用统一的方程表示：

P＝P₁+(P₂-P₁)H₁+(P₃-P₂)H₂

在Ω₁内，H₁＝H₂＝0；在Ω₂内，H₁＝1，H₂＝0；在Ω₃内，H₁＝H₂＝1。

交界面的位置确定方法，引入标记距离函数和指示函数H，或H＝0.5表示界面，重新初始化方程中包含界面的法线方向，由于指示函数H本身的局限性，不适宜用于求解远离H＝0.5处的界面法线方向，选取一个临界值δ，当时，利用计算界面法线方向，否则用H计算法线方向，具体过程表述如下：

(1)求解两个指示函数的输运方程，得到和H；

(2)对远离界面处的标记距离函数进行初始化；

(3)计算界面弥散区域各点的法线方向；

(4)对指示函数H进行初始化；

(5)根据H计算靠近界面处的距离标记函数

(6)进行物理特性插值，确定界面位置，进行下一步。

交界面的位置确定过程中选取的临界值δ一般可以取0.1Δ～0.2Δ，其中Δ表示网格尺寸，当时表示远离界面，表示靠近界面。

以步骤2中首浪产生数值模结果作为前提进行连贯的数值分析，用沿水深积分的Navier-Stokes方程即浅水方程作为控制方程求解涌浪的传播过程，通过步骤3的计算，可以得到各时刻各测点的水深数值，对涌浪的运动轨迹和泛滥区域进行数值预报，评价涌浪灾害。

典型滑坡体截面产生的初始浪高计算方法如下：在滑坡体进入水体之后的任意时刻，取指示函数H值分布图，取两点，分别为点A和点B，设H_A＜0.5，点A位于空气中，H_B＞0.5，点B位于水体中，A点高程h_A，B点高程h_B，根据线性插值远离，容易得出，H＝0.5的点的高程为初始浪高h＝h_O-h_Still，其中h_Still表示静水位。

本发明采用上述技术方案具有如下有益效果：

本发明针对库岸边坡失稳到涌浪产生与传播全过程进行数值模拟采用，将库岸边坡从失稳、滑动到涌浪产生再到涌浪传播看作是一个完整连续的过程，采用Navier-Stokes方程和改进的守恒式Level Set方法计算出首浪的各参数，在此基础上，采用浅水方程计算涌浪的传播，对涌浪运动轨迹和泛滥区域进行数值预报，进行涌浪灾害评价。

附图说明

图1是本发明三维有限元网格图，

图2是本发明典型断面密度分布图，

图3是本发明三相流指示函数H₁、H₂示意图，

图4是本发明首浪高度计算示意图，

图5是本发明河道截面处涌浪水位随时间变化过程(靠近滑坡体入水处)，

图6是本发明河道截面处涌浪水位随时间变化过程(靠近大坝)。

具体实施方式

下面结合附图和具体实施方式，进一步阐明本发明。应理解下述具体实施方式仅用于说明本发明而不用于限制本发明的范围。

本发明采用原始的不可压缩粘性流体Navier-Stokes方程描述水体、空气、滑坡体的相互作用。采用改进的守恒式Level Set方法，准确捕捉流体自由表面，真实模拟涌浪产生过程。在此基础上，采用沿水深积分的流体控制方程模拟涌浪传播过程，对涌浪运动轨迹和泛滥区域进行数值预报，评价涌浪灾害。

一种滑坡涌浪计算方法，包括如下步骤：

步骤1：根据已有地形、水位和滑坡体形状位置资料，建立三维有限元模型；

步骤2：结合Navier-Stokes方程和改进的守恒式Level Set方法，计算流体各运动参数，捕捉自由面；

步骤3：采用沿水深积分的流体控制方程模拟涌浪传播过程，对涌浪运动轨迹和泛滥区域进行数值预报。

步骤1具体如下：

根据实际地形资料，利用前处理软件，如GID、AutoCAD、Abaqus等建立水体、空气、滑坡体三维有限元模型。

步骤2具体如下：

步骤(2-1)对土质或松散介质坡体，基于“流体化”的特点，以流体力学基本方程Navier-Stokes方程配以非牛顿流体或拟流体本构关系描述其物理特性和运动过程：

进一步地，由于将滑坡体视作似流体，滑坡体、水体和空气三相介质可用同一本构模型表示：

式中，τ_y0表示流体屈服应力，φ为内摩擦角，I_2D是应变张量D的第二不变量，μ为黏性系数，η为表示流体特性的指数。

步骤(2-2)在利用Navier-Stokes方程计算出的流场信息的基础上，采用改进的守恒式双重Level Set方法捕捉滑坡体、库水和空气三者的交界面，反映三者间强烈的流固耦合作用，真实模拟滑坡体入水产生涌浪的过程。改进的守恒式双重Level Set方法包含两个部分：指示函数的输运，指示函数的重新初始化。

为了确定交界面的位置，引入标记距离函数和指示函数H。或H＝0.5表示界面。重新初始化方程中包含界面的法线方向，由于指示函数H本身的局限性，不适宜用于求解远离H＝0.5处的界面法线方向。选取一个临界值δ，当时，利用计算界面法线方向，否则用H计算法线方向。

步骤(2-2)的具体过程表述如下：

(1)在步骤(2-1)求解得到的速度场的基础上，求解两个指示函数的输运方程，得到和H；

(2)对远离界面处的标记距离函数进行初始化；

(3)计算界面弥散区域各点的法线方向；

(4)对指示函数H进行初始化；

(5)根据H计算靠近界面处的距离标记函数

(6)进行物理特性插值，确定界面位置，进行下一步。

进一步地，上述过程中选取的临界值δ一般可以取0.1Δ～0.2Δ，其中Δ表示网格尺寸，当时表示远离界面，表示靠近界面。

交替计算步骤(2-1)、(2-2)，可以得到三相流运动参数、捕捉界面位置，观察涌浪产生过程。选取初始涌浪作为涌浪传播计算的初始数值，用以计算后续的涌浪传播。

进一步地，在三相流中，各点的物理参数(如密度，黏性等)可以用指示函数H表示。引入指示函数H₁、H₂三相介质各自区域分别记作Ω₁、Ω₂、Ω₃，界面分别用Γ₁₂、Γ₂₃、Γ₁₃表示。P表示物理参数，则各区域参数可以用统一的方程表示：

P＝P₁+(P₂-P₁)H₁+(P₃-P₂)H₂

在Ω₁内，H₁＝H₂＝0；在Ω₂内，H₁＝1,H₂＝0；在Ω₃内，H₁＝H₂＝1。

步骤3具体如下：

以步骤2中首浪产生数值模结果作为前提进行连贯的数值分析，用沿水深积分的Navier-Stokes方程即浅水方程作为控制方程求解涌浪的传播过程。通过步骤3的计算，可以得到各时刻各测点的水深数值，对涌浪的运动轨迹和泛滥区域进行数值预报，评价涌浪灾害。

进一步地，典型滑坡体截面产生的初始浪高计算方法如下：

在滑坡体进入水体之后的任意时刻，取指示函数H值分布图。此处以二维示意图为例，三维情况同理。取两点，分别为点A，点B，不妨设H_A＜0.5(空气)，H_B＞0.5(水体)，A点高程h_A，B点高程h_B。根据线性插值远离，容易得出，H＝0.5的点的高程为初始浪高h＝h_O-h_Still，其中h_Still表示静水位。

下面举本发明的一个具体实施例：

某滑体干密度为2100kg/m³、饱和密度为2200kg/m³、水的密度为1000kg/m³、空气的密度为1kg/m³。

步骤1、根据已有地形、水位和滑坡体形状位置资料，做出合理的简化，建立三维有限元模型。

模型中地形表面坐标是从DEM数据获得，所以坐标系为大地坐标系。z轴为竖直向，向上为正。覆盖层滑坡涌浪模型结点41475个，单元207081个，其中滑坡体单元17261个。三维网格图如图1。边界条件有库岸边界和截取的库水边界。假设截取的库水为反射边界，即涌浪传到此处反射向回传播；库岸边界为反射边界，即涌浪传到岸坡上发生反射，然后继续向对岸传播并与入射波相互迭加和影响。涌浪传播河道视为天然河道，参考相关文献，粗糙系数n取0.05。

步骤2、结合Navier-Stokes方程和改进的守恒式Level Set方法，计算流体各运动参数，捕捉自由面。

(2-1)对土质或松散介质坡体，基于“流体化”的特点，以流体力学基本方程Navier-Stokes方程配以非牛顿流体或拟流体本构关系描述其物理特性和运动过程：

进一步地，由于将滑坡体视作似流体，三相介质可用同一本构模型表示：

式中，τ_y0表示流体屈服应力，φ为内摩擦角，I_2D是应变张量D的第二不变量，μ为黏性系数，η为表示流体特性的参量。

(2-2)在利用Navier-Stokes方程计算出的流场信息的基础上，采用改进的守恒式双重Level Set方法捕捉滑坡体、库水和空气三者的交界面，反映三者间强烈的流固耦合作用，真实模拟滑坡体入水产生涌浪的过程，典型断面密度分布图如图2。改进的守恒式双重Level Set方法包含两个部分：指示函数的输运，指示函数的重新初始化。

为了确定交界面的位置，引入标记距离函数和指示函数H。或H＝0.5表示界面。重新初始化方程中包含界面的法线方向，由于指示函数H本身的局限性，不适宜用于求解远离H＝0.5处的界面法线方向。选取一个临界值δ，当时，利用计算界面法线方向，否则用H计算法线方向。步骤(2-2)的具体过程表述如下：

(1)在步骤(2-1)求解得到的速度场的基础上，求解两个指示函数的输运方程，得到和H；

(2)对远离界面处的标记距离函数进行初始化；

(3)计算界面弥散区域各点的法线方向；

(4)对指示函数H进行初始化；

(5)根据H计算靠近界面处的距离标记函数

(6)进行物理特性插值，确定界面位置，进行下一步。

进一步地，上述过程中选取的临界值δ一般可以取0.1Δ～0.2Δ，其中Δ表示网格尺寸，当时表示远离界面，表示靠近界面。

交替计算步骤(2-1)、(2-2)，可以得到三相流运动参数、捕捉界面位置，观察涌浪产生过程，选取初始涌浪作为涌浪传播计算的初始数值，用以计算后续的涌浪传播。

进一步地，在三相流中，各点的物理参数(如密度，黏性等)可以用指示函数H表示。引入指示函数H₁,H₂，三相介质各自区域分别记作Ω₁,Ω₂,Ω₃，界面分别用Γ₁₂,Γ₂₃,Γ₁₃表示，如图3。P表示物理参数，则各区域参数可以用统一的方程表示：

P＝P₁+(P₂-P₁)H₁+(P₃-P₂)H₂

在Ω₁内，H₁＝H₂＝0；在Ω₂内，H₁＝1,H₂＝0；在Ω₃内，H₁＝H₂＝1。

步骤3、采用沿水深积分的流体控制方程模拟涌浪传播过程。

以步骤2中首浪产生数值模结果作为前提进行连贯的数值分析，用沿水深积分的Navier-Stokes方程即浅水方程作为控制方程求解涌浪的传播过程。通过步骤3的计算，可以得到各时刻各测点的水深数值，对涌浪的运动轨迹和泛滥区域进行数值预报，评价涌浪灾害。

进一步地，典型滑坡体截面产生的初始浪高计算方法如下：

在滑坡体进入水体之后的任意时刻，取指示函数H值分布图。此处以二维示意图为例，三维情况同理。在靠近界面处取两点，分别为点A，点B，不妨设H_A＜0.5(空气)，H_B＞0.5(水体)，A点高程h_A，B点高程h_B。根据线性插值远离，容易得出，H＝0.5的点的高程为初始浪高h＝h_O-h_Still，其中h_Still表示静水位。

取t＝22.13s时刻的涌浪高度作为首浪高度，如图4，计算结果见表1。

表1 典型断面涌浪计算参数及结果

在靠近滑坡体入水处和靠近大坝处选取两个截面，分别布置3个测点。两截面出水深随时间分布图分别见图5、图6。T＝700s时，涌浪传播到坝址附近并形成最大浪高为2.7m。故坝址处受涌浪影响最高水位为382.7m左右，比坝顶高程(384m)约低1.3m。计算得出的滑坡涌浪高度低于大坝高度，对大坝的冲击速度较低，在无外界不可抗力因素的影响下，滑坡涌浪对大坝存在一定的影响，但不会产生危害性的破坏。

综上所述，本发明提供了一种滑坡涌浪的计算方法，该方法将涌浪产生和传播作为一个统一连贯的过程研究，能表现出滑坡体、库水和空气三者强烈耦合作用。本发明不仅对库岸边坡失稳破坏及涌浪次生灾害的预测具有重大的理论意义和实用价值，而且对其它类似的滑坡、泥石流、堰塞湖等地质灾害的数值模拟亦具有一定的借鉴意义。

本发明方案所公开的技术手段不仅限于上述技术手段所公开的技术手段，还包括由以上技术特征任意组合所组成的技术方案。

以上述依据本发明的理想实施例为启示，通过上述的说明内容，相关工作人员完全可以在不偏离本项发明技术思想的范围内，进行多样的变更以及修改。本项发明的技术性范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其技术性范围。