CN113051634A - 一种库岸边坡滑坡涌浪的流固耦合数值计算方法 - Google Patents

Abstract

本发明公开了一种库岸边坡滑坡涌浪的流固耦合计算方法，本发明将基于欧拉框架的计算流体动力学方法与基于拉格朗日框架的离散元法耦合，采用守恒式Level Set方法捕捉流体自由表面，采用浸入边界法实现两种框架间物理量的转换，准确捕捉移动的物面边界，计算流固耦合作用力，避免由于使用经验公式计算流固耦合作用力造成的计算精度低的问题。本发明对流体域采用有限元建模，对于滑坡体采用离散元建模，通过迭代计算，实现滑坡体与流体间的真实耦合，尤其适用模拟大块岩石或岩石颗粒组成的滑坡体下滑造成的滑坡涌浪问题，在实际工程中有指导意义。

Description

一种库岸边坡滑坡涌浪的流固耦合数值计算方法

技术领域

本发明涉及水利工程灾害预测领域，具体涉及一种库岸边坡滑坡涌浪的流固耦合数值计算方法。

背景技术：

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

滑坡体运动的问题可以采用非连续介质力学方法进行模拟，其中离散元法由于其成熟性和显式求解的内在优势，得到了广泛应用。但如需模拟涌浪的产生过程，由于滑坡体、水体和空气间存在强烈的耦合作用，数值模拟的难度将大幅度增加。为了反映三相介质间的耦合作用，需要妥善描述滑坡体和流体间的界面以及水体和空气间的自由面，其中自由面的描述问题已经相对成熟，包括VOF法、Level Set法等。然而，对于移动的流固耦合界面的描述仍具有较大的挑战性，目前的方法依赖于过多的简化或假设。如多相流模型常用于模拟滑坡涌浪问题，但只适用于模拟含水率较高的滑坡体，然而，实际工程问题中，滑坡体常由大块岩石或岩石颗粒组成，不适宜用多相流模型进行模拟。因此，有学者提出采用不同方法模拟流体和滑坡体，如采用基于有限元框架的计算流体动力学模拟水和空气，利用离散元模拟滑坡体，但是，对于移动的流固耦合界面以及流固耦合作用力的计算则采用了较为简化的经验公式进行计算。此类经验公式只适用于圆形颗粒，将耦合作用力简化至颗粒形心处，使得颗粒的大小对流体不会造成影响，即流场不能准确感知颗粒的存在，因此，基于经验公式的流固耦合数值计算方法不能高精度的反映滑坡涌浪产生过程中的流固耦合作用。

尽管国内外学者对此已做了大量的研究工作，取得了丰富的研究成果，但滑坡体和流体间强烈耦合作用的准确数值模拟方法仍有待发展。综合目前的研究，本发明提供一种计算流体动力学与离散元耦合的数值模拟方法，采用浸入边界法实现两种框架间物理量的转换，解决滑坡体和流体间移动边界的问题，计算流固耦合作用力，避免了基于经验公式流固耦合方法计算精度低的问题，可用于模拟由大块岩石或岩石颗粒组成的滑坡体产生滑坡涌浪的过程。

发明内容

本发明的目的在于提供一种库岸边坡滑坡涌浪的流固耦合数值计算方法，以解决现有技术中导致的流固耦合数值计算方法不能高精度的反映滑坡涌浪产生过程中的流固耦合作用的缺陷。

一种库岸边坡滑坡涌浪动态自由面的计算方法，所述方法包括如下步骤：

获取预构建的边坡模型参数；

根据模型参数，计算得到流体的压力场以及流体的速度场；

根据流体的压力场计算滑坡体的运动速度；

对流体的压力场、流体的速度场以及滑坡体的运动速度进行收敛性判断；

满足收敛性后，根据速度场计算自由面位置；

根据自由面位置计算求解下一时刻自由面位置的流体域内各点物理参数。

进一步的，获取预构建的边坡模型参数的方法包括如下步骤：

根据已有地形、水位和滑坡体形状位置进行建模；

对流体域采用有限元建模；

对滑坡体采用离散元建模；

分别获取构建的流体域模型和滑坡体模型的模型参数。

进一步的，所述流体的压力场以及速度场的计算方法包括如下步骤：

在动量方程中引入附加体力项并计算其值；

根据附加体力项计算得到压力场增量；

通过压力场增量进行迭代更新得到压力场；

忽略动量方程中的附加体力项，计算得到中间速度场；

根据附加体力项对中间速度场进行校正得到速度场。

进一步的，在动量方程中引入附加体力项并计算其值的方法包括如下公式及步骤：

结合模型参数构建引入附加体力项的动量方程：

式中，u表示流体速度场；t表示时间；p表示流体压力场；ρ表示流体密度；f_b表示体力项；τ表示流体粘性应力张量；f表示附加体力项；

确定滑坡体外表面的速度边界条件：

Vⁿ⁺¹＝Uⁿ⁺¹；

式中，V表示物面边界点的固体速度；U表示由流体网格点插值到物面边界点上的速度；

根据速度边界条件确定附加体力项的计算公式：

fⁿ⁺¹Δt＝D(Fⁿ⁺¹Δt)＝D[Vⁿ⁺¹-I(uⁿ+Kⁿ+Rⁿ⁺¹)]；

式中：I表示插值函数，将流体网格点上的物理量插值到物面边界点上；D表示分配函数，将物面边界点上的物理量分配到流体网格点上；

Fⁿ⁺¹Δt＝Vⁿ⁺¹-I(uⁿ+Kⁿ+Rⁿ⁺¹)；

其中，Kⁿ由下式得到：

Rⁿ⁺¹由下式得到：

迭代求解体力项的计算公式：

令f^n+1,k,0＝f^n+1,k-1，i＝1，根据计算式：

f^n+1,k,i△t＝f^n+1,k,i-1△t+D[Vⁿ⁺¹-I(uⁿ+Kⁿ+Rⁿ⁺¹+f^n+1,k,i-1△t)]；

在此基础上，计算

进行收敛性判断；f^n+1,k-1表示第n+1个时间步的第k-1个外迭代步；f^n+1,k,0表示第n+1个时间步的第k-1个外迭代步里的内迭代步的初始值

其中||·||表示任一范数，如果范数小于容差ε，则结束，否则令i＝i+1，迭代计算直到满足收敛条件，最终得到：

f^n+1,k＝f^n+1,k,i。

进一步的，通过压力场增量进行迭代更新得到压力场的方法包括如下步骤：

pⁿ⁺¹＝pⁿ+△p^n,k；

式中，Δp^n,k为压力场增量根据

计算得到。

进一步的，对中间速度场进行校正得到速度场的方法为：

uⁿ⁺¹＝u^*+R^n+1,k；

式中，中间速度场u^*的计算式为u^*＝u^*,K+f^n+1,kΔt，u^*,K＝uⁿ+Kⁿ；u^*,k表示忽略了附加体力项的中间速度场；f^n+1,k△t表示附加体力和时间的乘积，具有速度的量纲。

进一步的，所述滑坡体的相互作用力的计算方法包括如下步骤：

根据牛顿第二定律计算滑坡体的运动，控制方程如下：

式中，m表示离散滑坡体质量；

表示平动加速度；

上式中，F表示作用在浸入边界点上的力，F_e表示施加在离散滑坡体上的外力；F_c表示离散滑坡体间的相互作用；F_f表示流体作用在离散滑坡体上的力；c表示离散滑坡体的阻尼系数，

表示平动速度；

流体作用于离散滑坡体上的力F_f的计算公式：

式中，Ω′表示整个计算域；包括流体和离散滑坡体内部区域；Ω″表示离散滑坡体内部区域；f_b表示体力项；n表示计算时间步，k表示迭代时间步。

进一步的，对流体的压力场、流体的速度场以及滑坡体的运动速度进行收敛性判断的方法包括如下步骤：

根据下式：

其中

可取为u，p以及δ，如果计算结果小于

则满足收敛性检查，否则重复计算u，p以及δ，直到满足收敛性判断条件，

表示容差。

进一步的，所述自由面位置的运算方法包括如下步骤：

计算流体域指示函数的值：

确定对流方程：

初始化对流方程：

式中，τ表示虚拟时间步；n表示自由面法向；H表示自由面。

进一步的，根据涌浪自由面获取下一时刻的参数的方法包括如下步骤：

寻找自由面位置：

H＝0.5；

根据自由面位置，计算流体域内各点物理参数；

物理参数MP可以用指示函数H表示，空气域记作Ω₁，物理参数记为MP₁，水体域记作Ω₂，物理参数记为MP₂；

当待求点位于Ω₁时，H＝0，当待求点位于Ω₂时，H＝1；

流体计算域内各点的物理参数MP计算式如下：

MP＝MP₁+(MP₂-MP₁)H。

本发明的优点在于：该种库岸边坡滑坡涌浪动态自由面的计算方法：

(1)将基于欧拉框架的计算流体动力学方法与基于拉格朗日框架的离散元法耦合，采用浸入边界法实现两种框架间物理量的转换，准确捕捉移动的物面边界，计算流固耦合作用力，避免由于使用经验公式计算流固耦合作用力造成的计算精度低的问题；

(2)对流体域采用有限元建模，对于滑坡体采用离散元建模，尤其适用模拟滑坡体下滑造成的滑坡涌浪问题，在实际工程中有指导意义。

附图说明

图1是实施例一颗粒沉降初始时刻示意图；

图2是实施例一基于经验公式法的流固耦合法在t＝1s时的计算结果；

图3是实施例一基于本发明提出流固耦合法在t＝0.5s时的计算结果；

图4是实施例一基于本发明提出流固耦合法在t＝1.51s时的计算结果；

图5是实施例一基于本发明提出流固耦合法在t＝15.04s时的计算结果；

图6是实施例二滑坡涌浪初始时刻示意图；

图7是实施例二t＝0.76s时滑坡涌浪示意图；

图8是实施例二t＝5.95s时滑坡涌浪示意图；

图9是实施例二t＝7.68s时滑坡涌浪示意图；

图10是实施例二t＝11.14s时滑坡涌浪示意图。

具体实施方式

为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施方式，进一步阐述本发明。

如图1至图10所示，一种库岸边坡滑坡涌浪动态自由面的计算方法，其特征在于，所述方法包括如下步骤：

获取预构建的边坡模型参数；

根据模型参数，计算得到流体的压力场以及流体的速度场；

根据流体的压力场计算滑坡体的运动速度；

对流体的压力场、流体的速度场以及滑坡体的运动速度进行收敛性判断；

满足收敛性后，根据速度场计算自由面位置；

根据自由面位置计算求解下一时刻自由面位置的流体域内各点物理参数。

具体步骤包括：

步骤1：

根据已有地形、水位和滑坡体形状位置资料，建立数值计算模型，对流体域采用有限元建模，对于滑坡体采用离散元建模。

步骤2：

以不可压缩流体Navier-Stokes方程作为控制方程，为了表示滑坡体和流体间的移动界面以及考虑滑坡体对流体的作用力，结合浸入边界法的思想，在动量方程右端引入附加体力项f，最终，流体的控制方程形式如下：

连续性方程：

动量方程：

式中，u表示流场速度，t表示时间，p表示压力，ρ表示密度，f_b表示体力项，τ表示粘性应力张量，μ表示动力粘性系数。

根据滑坡体外表面的速度边界条件，Vⁿ⁺¹＝Uⁿ⁺¹，其中V表示物面边界点的固体速度，U表示由流体网格点插值到物面边界点上的速度，可以得出附加体力项的计算公式

fⁿ⁺¹Δt＝D(Fⁿ⁺¹Δt)＝D[Vⁿ⁺¹-I(uⁿ+Kⁿ+Rⁿ⁺¹)]

式中，I表示插值函数，将流体网格点上的物理量插值到物面边界点上；D表示分配函数，将物面边界点上的物理量分配到流体网格点上，本发明中采用Peskin(2002)给出的插值、分配函数；Kⁿ、

的表达式如下：

步骤2.1忽略压力及附加体力项，计算中间速度场u^*,K

中间速度场u^*,K的计算式为u^*,K＝uⁿ+Kⁿ。

步骤2.2计算中间压力场Δp^n,k及附加体力项f^n+1,k

令p^n+1,0＝pⁿ，f^n+1,0＝0，k＝1。

步骤2.2.1根据附加体力项的计算公式计算f^n+1,k

由于浸入边界点与网格点不重合，因此通过单次计算附加体力项的计算公式，并不能满足速度边界条件，所以需要迭代求解体力项的计算公式。令f^n+1,k,0＝f^n+1,k-1，i＝1，根据计算式f^n+1,k,i△t＝f^n+1,k,i-1△t+D[Vⁿ⁺¹-I(uⁿ+Kⁿ+Rⁿ⁺¹+f^n+1,k,i-1△t)]计算f^n+1,k,i。在此基础上，计算

进行收敛性判断，其中||·||表示任一范数，如果范数小于容差ε，则结束步骤2.2.1，否则令i＝i+1，迭代计算直到满足收敛条件。最终，f^n+1,k＝f^n+1,k,i。

求f的时候，是需要迭代的，把这个迭代叫做内迭代，i表示这个迭代的第i步；流体的速度u，压力p以及f是耦合的，需要迭代求解，这个迭代叫做外迭代，k表示这个迭代的第k步；每个外迭代里面都需要进行内迭代求f；f^n+1,k-1表示第n+1个时间步的第k-1个外迭代步；

f^n+1,k,0表示第n+1个时间步的第k-1个外迭代步里的内迭代步的初始值；

步骤2.2.2更新中间速度场u^*

更新中间速度场u^*的计算式为u^*＝u^*,K+f^n+1,kΔt。

步骤2.2.3根据连续性方程压力场增量Δp^n,k

根据

计算压力场增量Δp^n,k。

步骤2.2.4收敛性检查

计算||I(uⁿ+Kⁿ+R^n+1,k)-I(uⁿ+Kⁿ+R^n+1,k-1)||，如果小于容差ε，则满足收敛性检查，结束2.2步的迭代计算，否则，令k＝k+1，重复步骤2.2.1-步骤2.2.4，直到满足收敛条件。

步骤2.3更新压力场pⁿ⁺¹，校正速度场uⁿ⁺¹

pⁿ⁺¹＝pⁿ+△p^n,k，

uⁿ⁺¹＝u^*+R^n+1,k。

步骤3：

根据牛顿第二定律计算滑坡体的运动，控制方程如下

式中，m表示质量，

表示平动加速度，

F_e表示施加在离散滑坡体上的外力，F_c表示离散滑坡体间的相互作用，F_f表示流体作用在离散滑坡体上的力，c表示粘性系数，

表示平动速度。

步骤3.1计算离散滑坡体间的相互作用力F_c

根据赵兰浩等(ZL201610218708.3)提出的距离势函数离散单元法计算滑坡体间的相互作用力。

步骤3.2计算流体作用在离散滑坡体上的力F_f

流体作用在离散滑坡体上的力的计算式为

式中，F表示作用在浸入边界点上的力，Ω″表示离散滑坡体内部区域。

步骤4：

计算

其中

可取为u，p以及δ，如果小于

则满足收敛性检查，执行步骤5，

表示容差。否则重复步骤2-步骤4，直到满足收敛性判断条件。

步骤5：

步骤5.1根据CLS方法计算流体域指示函数H的值，包括

对流方程：

重新初始化方程：

式中，τ表示虚拟时间步，n表示自由面法向。

步骤5.2寻找自由面位置

H＝0.5表示自由面

步骤5.3计算流体域内各点物理参数

流体计算域内各点的物理参数MP(如密度，黏性等)可以用指示函数H表示，空气域记作Ω₁，物理参数记为MP₁，水体域记作Ω₂，物理参数记为MP₂，当待求点位于Ω₁时，H＝0，当待求点位于Ω₂时，H＝1。流体计算域内各点的物理参数MP计算式如下

MP＝MP₁+(MP₂-MP₁)H。

步骤6：

重复步骤2至步骤5，直至计算完所有时间步。

以下提供两个实施例进一步说明：

实施例一：

多颗粒沉降数值计算：初始情况如图1，计算域大小为10m×10m，64个颗粒可分为两类，直径分别为0.550m、0.275m。流体密度与颗粒密度分别为1000kg/m³、2600kg/m³。颗粒的法向和切向刚度均为4.8×10⁸Pa。颗粒由静止状态开始，在重力作用下沉降。采用基于经验公式的流固耦合方法计算得到的t＝1s的结果如图2。采用本发明中流固耦合方法，将流体域划分为480×480个规则四边形有限元单元，在颗粒表面分别布置584、1551个浸入边界点，计算得到的沉降过程如图3-图5，时间分别为t＝0.5s、t＝1.51s、t＝15.04s。

可以观察到，采用基于经验公式的流固耦合方法，得到的流场不能感知到颗粒的存在，颗粒的大小对流场没有影响。采用本发明提供的方法，准确捕捉移动的物面边界，提高流场及颗粒运动的计算精度。

实施例二:

初始时刻，滑坡体静止，位置如图6所示，滑坡体和水体的密度分别为2640kg/m³和1000kg/m³，边坡的内摩擦角为43°，动床摩擦角为24°。采用本发明提供的方法，将流体域划分为753664个不规则三角形有限元单元，将滑坡体划分为2048个三角形离散元单元。

图7-图10为四个时刻(t＝0.76s、t＝5.95s、t＝7.68s、t＝11.14s)滑坡涌浪产生过程示意图。

可以观察到，采用本发明提供的方法，可以清晰模拟涌浪产生过程。

由技术常识可知，本发明可以通过其它的不脱离其精神实质或必要特征的实施方案来实现。因此，上述公开的实施方案，就各方面而言，都只是举例说明，并不是仅有的。所有在本发明范围内或在等同于本发明的范围内的改变均被本发明包含。

Claims (10)

1.一种库岸边坡滑坡涌浪动态自由面的计算方法，其特征在于，所述方法包括如下步骤：

获取预构建的边坡模型参数；

根据模型参数，计算得到流体的压力场以及流体的速度场；

根据流体的压力场计算滑坡体的运动速度；

对流体的压力场、流体的速度场以及滑坡体的运动速度进行收敛性判断；

满足收敛性后，根据速度场计算自由面位置；

根据自由面位置计算求解下一时刻自由面位置的流体域内各点物理参数。

2.根据权利要求1所述的一种库岸边坡滑坡涌浪动态自由面的计算方法，其特征在于：获取预构建的边坡模型参数的方法包括如下步骤：

根据已有地形、水位和滑坡体形状位置进行建模；

对流体域采用有限元建模；

对滑坡体采用离散元建模；

分别获取构建的流体域模型和滑坡体模型的模型参数。

3.根据权利要求1所述的一种库岸边坡滑坡涌浪动态自由面的计算方法，其特征在于：所述流体的压力场以及速度场的计算方法包括如下步骤：

在动量方程中引入附加体力项并计算其值；

根据附加体力项计算得到压力场增量；

通过压力场增量进行迭代更新得到压力场；

忽略动量方程中的附加体力项，计算得到中间速度场；

根据附加体力项对中间速度场进行校正得到速度场。

4.根据权利要求3所述的一种库岸边坡滑坡涌浪动态自由面的计算方法，其特征在于：在动量方程中引入附加体力项并计算其值的方法包括如下公式及步骤：

结合模型参数构建引入附加体力项的动量方程：

式中，u表示流体速度场；t表示时间；p表示流体压力场；ρ表示流体密度；f_b表示体力项；τ表示流体粘性应力张量；f表示附加体力项；

确定滑坡体外表面的速度边界条件：

Vⁿ⁺¹＝Uⁿ⁺¹；

式中，V表示物面边界点的固体速度；U表示由流体网格点插值到物面边界点上的速度；

根据速度边界条件确定附加体力项的计算公式：

fⁿ⁺¹Δt＝D(Fⁿ⁺¹Δt)＝D[Vⁿ⁺¹-I(uⁿ+Kⁿ+Rⁿ⁺¹)]；

式中：I表示插值函数，将流体网格点上的物理量插值到物面边界点上；D表示分配函数，将物面边界点上的物理量分配到流体网格点上；

Fⁿ⁺¹Δt＝Vⁿ⁺¹-I(uⁿ+Kⁿ+Rⁿ⁺¹)；

其中，Kⁿ由下式得到：

Rⁿ⁺¹由下式得到：

迭代求解体力项的计算公式：

令f^n+1,k,0＝f^n+1,k-1，i＝1，根据计算式：

f^n+1,k,i△t＝f^n+1,k,i-1△t+D[Vⁿ⁺¹-I(uⁿ+Kⁿ+Rⁿ⁺¹+f^n+1,k,i-1△t)]；

在此基础上，计算

进行收敛性判断；f^n+1,k-1表示第n+1个时间步的第k-1个外迭代步，f^n+1,k,0表示第n+1个时间步的第k-1个外迭代步里的内迭代步的初始值。

其中||·||表示任一范数，如果范数小于容差ε，则结束，否则令i＝i+1，迭代计算直到满足收敛条件，最终得到：

f^n+1,k＝f^n+1,k,i。

5.根据权利要求4所述的一种库岸边坡滑坡涌浪动态自由面的计算方法，其特征在于：通过压力场增量进行迭代更新得到压力场的方法包括如下步骤：

pⁿ⁺¹＝pⁿ+△pⁿ,^k；

式中，Δp^n,k为压力场增量根据

计算得到。

6.根据权利要求5所述的一种库岸边坡滑坡涌浪动态自由面的计算方法，其特征在于：对中间速度场进行校正得到速度场的方法为：

uⁿ⁺¹＝u^*+R^n+1,k；

式中，中间速度场u^*的计算式为u^*＝u^*,K+f^n+1,kΔt，u^*,K＝uⁿ+Kⁿ；u^*,k表示忽略了附加体力项的中间速度场；f^n+1,k△t表示附加体力和时间的乘积，具有速度的量纲。

7.根据权利要求6所述的一种库岸边坡滑坡涌浪动态自由面的计算方法，其特征在于：所述滑坡体的相互作用力的计算方法包括如下步骤：

根据牛顿第二定律计算滑坡体的运动，控制方程如下：

式中，m表示离散滑坡体质量；

表示平动加速度；

上式中，F表示作用在浸入边界点上的力，F_e表示施加在离散滑坡体上的外力；F_c表示离散滑坡体间的相互作用；F_f表示流体作用在离散滑坡体上的力；c表示离散滑坡体的阻尼系数，

表示平动速度；

流体作用于离散滑坡体上的力F_f的计算公式：

式中，Ω′表示整个计算域；包括流体和离散滑坡体内部区域；Ω″表示离散滑坡体内部区域；f_b表示体力项；n表示计算时间步，k表示迭代时间步。

8.根据权利要求7所述的一种库岸边坡滑坡涌浪动态自由面的计算方法，其特征在于：对流体的压力场、流体的速度场以及滑坡体的运动速度进行收敛性判断的方法包括如下步骤：

根据下式：

其中

可取为u，p以及δ，如果计算结果小于

则满足收敛性检查，否则重复计算u，p以及δ，直到满足收敛性判断条件，

表示容差。

9.根据权利要求8所述的一种库岸边坡滑坡涌浪动态自由面的计算方法，其特征在于：所述自由面位置的运算方法包括如下步骤：

计算流体域指示函数的值：

确定对流方程：

初始化对流方程：

式中，τ表示虚拟时间步；n表示自由面法向；H表示自由面。

10.根据权利要求9所述的一种库岸边坡滑坡涌浪动态自由面的计算方法，其特征在于：根据自由面位置计算求解下一时刻自由面位置的流体域内各点物理参数的方法包括如下步骤：

寻找自由面位置：

H＝0.5；

根据自由面位置，计算流体域内各点物理参数；

物理参数MP可以用指示函数H表示，空气域记作Ω₁，物理参数记为MP₁，水体域记作Ω₂，物理参数记为MP₂；

当待求点位于Ω₁时，H＝0，当待求点位于Ω₂时，H＝1；

流体计算域内各点的物理参数MP计算式如下：

MP＝MP₁+(MP₂-MP₁)H。

Images