一种基于FDM-DDA的非饱和瞬态流固耦合计算方法

摘要

本发明涉及数据处理领域，具体涉及一种基于FDM‑DDA的非饱和瞬态流固耦合计算方法，包括：采用matlab生成裂隙以建立岩土体模型；对所述岩土体模型中的岩体裂隙进行识别，形成流体网格；对流体网络中的岩体裂隙开度计算；进行DDA非饱和稳态渗流场的计算，利用FDM进行非饱和瞬态渗流场计算，得到稳态计算结果；将稳态计算结果作为流体计算模块写入DDA并赋予非饱和瞬态渗流场计算参数，求解各节点处的水头和体积含水量作为瞬态计算结果；基于瞬态计算结果进行流固耦合分析。本发明增加了岩体裂隙识别的功能，利用FDM将richards方程写入DDA中，实现了DDA对岩土体非饱和瞬态渗流过程的模拟，能够用于降雨、库水位变化等多种工况中。

说明书

技术领域

本发明涉及数据处理领域，尤其涉及一种基于FDM-DDA的非饱和瞬态流固耦合计算方法。

背景技术

在对边坡工程进行安全稳定性分析时，往往会考虑暴雨工况下边坡的状态。降雨会对岩土体的强度起到弱化作用，因此降雨往往是边坡失稳的重要诱发因素，也对边坡带来了安全隐患。

目前，针对边坡工程的岩土体的建模中，有限单元法被广泛应用，但是有限单元法基于连续性假设，不符合岩土体作为离散介质的本质。因此DDA作为一种具有大变形大位移特点的数值方法被广泛应用于岩土体的数值模拟中，以 DDA为平台的仿真对工程设计、施工稳定性评估起到了积极的意义。

但是其存在以下弊端：

1.DDA的岩土体数值仿真多集中于固体方面，即单纯的岩石与土体方面，忽略了渗流对岩土体的影响；

2.DDA在考虑水-岩土体相互作用时，即流固耦合模拟过程中，往往只考虑了稳态渗流。对于降雨、库水位变动等非稳定的瞬态过程，采用稳态渗流进行模拟并不能反应其变化过程对岩土体造成的影响。

针对以上问题，亟需一种适用于DDA的瞬态流固耦合计算方法，以实现 DDA对降雨、库水位变动等非稳定渗流岩土体的模拟。

发明内容

本发明的目的在于提供一种基于FDM-DDA的非饱和瞬态流固耦合计算方法，旨在利用DDA进行前处理并进行应力场计算。利用FDM进行非饱和瞬态渗流场计算，再将其结果导入DDA进行流固耦合分析，以实现用DDA对降雨、库水位变动等非稳定渗流岩土体的模拟。

为实现上述目的，本发明提供了一种基于FDM-DDA的非饱和瞬态流固耦合计算方法，包括采用matlab生成裂隙以建立岩土体模型；

对所述岩土体模型中的岩体裂隙进行识别，形成流体网格；

对流体网络中的岩体裂隙开度计算；

进行DDA非饱和稳态渗流场的计算，利用FDM进行非饱和瞬态渗流场计算，得到稳态计算结果；

将稳态计算结果作为流体计算模块写入DDA并赋予非饱和瞬态渗流场计算参数，求解各节点处的水头和体积含水量作为瞬态计算结果；

基于瞬态计算结果进行流固耦合分析。

其中，所述采用matlab生成裂隙的具体方式是：采用matlab编写Monte Carlo 算法，生成随机裂隙的两端点坐标，将其导入DDA的前处理dc程序以生成裂隙。

其中，所述对岩土体模型中的岩体裂隙进行识别，形成流体网格的具体步骤是：

获取相邻两点形成的线段并将端点顺序标号，并将标号存储到第一二维数据的第一列；

对所述线段进行遍历，寻找端点坐标的差值在预设范围内且不属于同一块体的两条线段作为短裂隙的两侧；

将短裂隙的两侧的线段序号存储在第一二维数据的第二列中；

在所述线段中基于平行线的方式确认长裂隙；

将长裂隙拆分形成短裂隙，并将短裂隙的线段序号存储到第二二维数组的前两列中；

将第一二维数组和第二二维数组中的相同点进行合并以形成流体网络。

其中，所述对流体网络中的岩体裂隙开度计算的具体步骤是：

获取裂隙两侧线段端点序号，一侧的序号为第一点、第二点，另一侧为第三点和第四点；

基于第一点、第二点计算到第三点的第一距离；

基于第一点、第二点计算到第四点的第二距离；

基于第一距离和第二距离计算等效水力开度。

其中，所述进行DDA非饱和稳态渗流场的计算的具体步骤是：

采用立方定律对单裂隙渗流进行计算；

采用Gauss–Seidel迭代求得每个节点的水头，再计算每个节点的体积含水量；

将相关程序作为稳态渗流模块写入DDA程序中。

其中，所述利用FDM进行非饱和瞬态渗流场计算的具体步骤是：

将岩体裂隙作为FDM渗流场计算所需的流体网格；

采用FDM将richards方程转化成符合DDA格式的体积含水量矩阵。

其中，所述基于瞬态计算结果进行流固耦合分析的具体步骤是：

将瞬态计算结果根据最小势能原理转化成水荷载矩阵；

将水载荷矩阵加入到DDA的荷载矩阵中进行DDA应力场计算，当块体位移小于规定值时，应力场计算达到平衡。

本发明的一种基于FDM-DDA的非饱和瞬态流固耦合计算方法，包括：采用matlab生成裂隙以建立岩土体模型；对所述岩土体模型中的岩体裂隙进行识别，形成流体网格；对流体网络中的岩体裂隙开度计算；进行DDA非饱和稳态渗流场的计算，利用FDM进行非饱和瞬态渗流场计算，得到稳态计算结果；将稳态计算结果作为流体计算模块写入DDA并赋予非饱和瞬态渗流场计算参数，求解各节点处的水头和体积含水量作为瞬态计算结果；基于瞬态计算结果进行流固耦合分析。本发明针对DDA无法对岩土体非饱和瞬态渗流进行模拟的缺点进行了改进。增加了岩体裂隙识别的功能，能对凹多边形、长裂隙等较复杂情况进行识别，对裂隙开度进行计算，以及对模型外边界进行区分；本发明利用 FDM将richards方程写入DDA中，实现了DDA对岩土体非饱和瞬态渗流过程的模拟，能够用于降雨、库水位变化等多种工况中；本发明能够考虑岩土体与水体的相互作用，即流固耦合计算，更加符合工程中降雨等流体对岩土工程影响的工程实际，进一步地为工程设计、安全评估提供了益处。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明的本发明DDA随机裂隙岩土体计算模型的结构图。

图2是本发明的本发明岩体裂隙识别得到二维数组的输出文件示意图；

图3是本发明的本发明流固耦合计算结果有效饱和度云图；图4是图1的剖面示意图。

图4是本发明的一种基于FDM-DDA的非饱和瞬态流固耦合计算方法的流程图。

图5是本发明的对所述岩土体模型中的岩体裂隙进行识别，形成流体网格的流程图。

图6是本发明的对流体网络中的岩体裂隙开度计算的流程图。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

请参阅图1～图6，本发明提供一种基于FDM-DDA的非饱和瞬态流固耦合计算方法，包括：

S101采用matlab生成裂隙以建立岩土体模型；

所述采用matlab生成裂隙的具体方式是：采用matlab编写Monte Carlo算法，生成随机裂隙的两端点坐标，将其导入DDA的前处理dc程序以生成裂隙。

根据统计学原理确定岩体裂隙中点、倾向、倾角、裂长符合的概型分布，利用MATLAB的随机数种子通过random命令生成符合指定概型的随机数。根据随机生成裂隙的直线点斜式求解裂隙两端点的坐标。将坐标数据输出，通过 AUTOCAD脚本文件.scr进行绘图。将绘图结果储存成.dxf格式文件，分线和点两种图层导入DDA实现岩土体模型建立。

S102对所述岩土体模型中的岩体裂隙进行识别，形成流体网格；

具体步骤是：

S201获取相邻两点形成的线段并将端点顺序标号，并将标号存储到第一二维数据的第一列；

DDA的块体由顶点坐标连线形成，从块体一点出发逆时针连线回到初始点。根据DDA数据的特点，将点i与下一个点i+1形成的线段用点i的序号表示，并存储在第一二维数组lem的第一列中。

S202对所述线段进行遍历，寻找端点坐标的差值在预设范围内且不属于同一块体的两条线段作为短裂隙的两侧；

对线段进行遍历，寻找端点坐标足够接近且不属于同一块体的两条线段作为岩体裂隙的两侧。将寻找到的线段序号存储在二维数组的第二列中。

S203将短裂隙的两侧的线段序号存储在第一二维数据的第二列中；

将寻找到的线段序号存储在二维数组的第二列中。

S204在所述线段中基于平行线的方式确认长裂隙；

在长裂隙中，DDA块体可能在裂隙的一侧没有块体顶点，因此上述方式不能识别出该类型的裂隙。除开这种形式的裂隙，其他的都是短裂隙。因为每个序号i只唯一对应由点i和i+1组成的线段，所以长裂隙只需要在没有存储到第一二维数组的节点中进行筛选。

S205将长裂隙拆分形成短裂隙，并将短裂隙的线段序号存储到第二二维数组的前两列中；

对于这一类型的节点，通过寻找平行线的方式来确定该裂隙的另一侧。找到平行线后，根据线段的包含关系来将长裂隙拆分成几条短裂隙组合成的裂隙。对于拆分形成的短裂隙通过向量的夹角余弦来判断是否还有点落在了短裂隙的线段上。若有，则对该裂隙进一步的拆分；若无，则将该线段的两端点序号存入第二二维数组lemm的前两列中。

S206将第一二维数组和第二二维数组中的相同点进行合并以形成流体网络。

最后将在以上两个二维数组中将相同点进行合并。分别在lem和lemm数组内，以及lem和lemm数组之间进行合并。当两点的坐标足够接近认为是同一个点，代表了经过同一个节点的两条不同线段。通过合并线段找到每一个节点与其周围节点的连接关系，并将合并的结果储存进二维数组leq的前两列中。

S103对流体网络中的岩体裂隙开度计算；

具体步骤是：

S301获取裂隙两侧线段端点序号，一侧的序号为第一点、第二点，另一侧为第三点和第四点；

S302基于第一点、第二点计算到第三点的第一距离；

根据上面的步骤中储存的岩体裂隙两侧线段序号，通过公式(1)计算裂隙一侧到另一侧其中一个端点的第一距离e

其中，下标1、2表示裂隙同一侧的两端点坐标(x1,y1),(x2,y2)，下标 3表示裂隙另一侧端点的坐标(x3,y3)。

S303基于第一点、第二点计算到第四点的第二距离；

同理可以计算端点1、2与另一侧端点4的第二距离e

S304基于第一距离和第二距离计算等效水力开度。

根据经验公式来计算等效水力开度。

式(2)中r是距离e

S104进行DDA非饱和稳态渗流场的计算，利用FDM进行非饱和瞬态渗流场计算，得到稳态计算结果；

所述进行DDA非饱和稳态渗流场的计算的具体步骤是：

401采用立方定律对单裂隙渗流进行计算；

采用立方定律对单裂隙渗流进行计算，按照节点的位置关系确定渗流方向并将单裂隙的流量叠加可以得到式(3)

式(3)中T

S402采用Gauss–Seidel迭代求得每个节点的水头，再计算每个节点的体积含水量；

采用Gauss–Seidel迭代求得每个节点的水头，再计算每个节点的体积含水量。

S

式(4)中P

S403将相关程序作为稳态渗流模块写入DDA程序中。

所述利用FDM进行非饱和瞬态渗流场计算的具体步骤是：

S501将岩体裂隙作为FDM渗流场计算所需的流体网格；

S502采用FDM将richards方程转化成符合DDA格式的体积含水量矩阵。

在二维状态下的richards方程如下所示：

式(5)中K(θ)是非饱和渗流系数，是与岩体的体积含水量θ相关的函数。将其改写成体积含水量的表达形式。

式(6)中

式(7)中θ上标表示时间差分，下标表示空间差分。对y方向也采用相同的方式计算。继续对式(7)中的偏导进行离散化。

由于有限差分在非饱和渗流中存在一定的滞后性，采用D(θ

对y方向也有类似的计算过程，为了简化表达令

式(11)中上标为j表示本时步的状态，均为已知量。j+1表示下一时步的状态，通过本时步的状态求解。Δt为瞬态渗流的时间步长，将式(11)整理成矩阵形式。

式(12)中A

A

S105将稳态计算结果作为流体计算模块写入DDA并赋予非饱和瞬态渗流场计算参数，求解各节点处的水头和体积含水量作为瞬态计算结果；

式(12)系数矩阵A是一个严格对角占优的稀疏矩阵，与DDA类似的采用 SOR迭代求解。要求解式(12)需要求得时步j对应的体积含水量。

式(13)中S

S106基于瞬态计算结果进行流固耦合分析。

具体步骤是：

S601将瞬态计算结果根据最小势能原理转化成水荷载矩阵；

S602将水载荷矩阵加入到DDA的荷载矩阵中进行DDA应力场计算，当块体位移小于规定值时，应力场计算达到平衡。

DDA的总体平衡方程如下：

式(14)中K是块体的刚度矩阵，D是变形矩阵，F是荷载矩阵。将步骤 S105得到的节点水头根据最小势能原理转化成水荷载矩阵，加入到DDA的荷载矩阵F中进行DDA应力场计算。当某一时步的块体位移小于规定值时，认为应力场计算达到平衡。更新步骤S103中裂隙开度并进行渗流模块计算，然后重复步骤S104和S105，以实现DDA流固耦合分析。

下面以二维问题为例对发明进一步进行说明。建立岩土体模型并导入DDA。设置两组裂隙，其倾向分别为100～128°、180～185°，倾角为73～84°、69～ 78°，裂隙长为3～5m、2～30m，由编写的matlab程序生成随机裂隙并导入 DDA中得到图1；经过步骤2的岩体裂隙识别得到的数组如图2所示，设置岩体的初始裂隙开度为0.0003m.当步骤S103计算得到的等效水力开度大于该值时进行更新。DDA非饱和瞬态渗流场计算。将模型的底部设置为0水头边界，模型顶部在降雨量小于渗透系数时设置为流量边界，即在式(3)中将对应点的设置成已知量。降雨量小于渗透系数时设置为水头边界，即在式(3)中将对应点的设置成已知量，将设置为1其余的设置为0。再设置体积含水量对应的初始条件，即在式(12)中将每个节点的设置为步骤S104稳态渗流场计算对应的结果，边界上的设置成对应的已知量，即可实现DDA非饱和瞬态渗流场计算。

以上所揭露的仅为本发明一种较佳实施例而已，当然不能以此来限定本发明之权利范围，本领域普通技术人员可以理解实现上述实施例的全部或部分流程，并依本发明权利要求所作的等同变化，仍属于发明所涵盖的范围。