基于钻孔数据的地层力学参数有限元建模方法及终端设备

1.本发明涉及有限元建模领域，特别是一种基于钻孔数据的地层力学参数有限元建模方法及终端设备。


背景技术：

2.随着计算机技术的进步与发展，有限元方法在土木工程设计、计算等方面发挥着越来越重要的作用。关于地下结构的有限元计算模型主要分为基于松弛荷载理论的结构力学模型和基于岩承理论的岩体力学模型两类，其中岩体力学模型在有限元建模时需要将地下结构与周围地层共同建立起来。在岩体力学模型中地层建模时能否准确反映地层现实特性具有重要意义。
3.现有的建模方法往往是在钻孔数据所得的地质图选取其中某一个或几个受力不利断面进行建模计算，具体在建模计算时则忽略地层的起伏特性，依据地层平均厚度将地层建成水平状，随后进行地下结构的设计。这种设计方法应用于地层起伏变化不大的地层尚可，但应用于地层起伏变化较大的地层中有两方面弊端，1)在横断面结构设计中无法考虑地层起伏变化所带来的偏压效果，使得计算结果不准确；2)在纵断面设计计算中，无法体现地层起伏变化所导致的纵向不均匀荷载。例如在计算盾构隧道顶推力时，若将地层考虑为水平状，则所得顶推力结果为定值，显然与实际情况不符。因此，在地层起伏变化较大情况下有限元建模计算时必须考虑地层沿纵向的不均匀分布。
4.现有的考虑地层起伏特性的有限元建模主要有两种方法，一种是在几何建模时将地层起伏形状建出，另一种则是考虑地层力学参数空间变异性的基于随机场理论的方法(例如cn113158315a)。其中第一种方法常会遇到某些地层形状出现尖端情况，此时在网格划分时极易遇到网格畸变等问题，易导致计算结果误差较大。而第二种方法虽可以表征地层分布的不均匀性，但具体实现过程较为复杂，且需要大量的计算才能得到统计结果的稳定值，计算量大。


技术实现要素：

5.本发明所要解决的技术问题是，针对现有技术不足，提供一种基于钻孔数据的地层力学参数有限元建模方法及终端设备，避免随机场方法庞大的计算量。避免从几何建模角度在有限元模型直接划分地层所带来的网格畸变问题。
6.为解决上述技术问题，本发明所采用的技术方案是：一种基于钻孔数据的地层力学参数有限元建模方法，包括以下步骤：
7.s1、建立无地层参数的有限元模型；采集工程现场钻孔数据；所述工程现场钻孔数据包括各钻孔位置及编号、各钻孔内各土层的埋深及力学参数；
8.s2、离散化所述现场钻孔数据，并对所述有限元模型赋予相应参数；
9.s3、根据有限元模型单元间的相关作用关系，计算非钻孔位置单元的参数；
10.s4、更新并分类非钻孔位置地层参数，并将更新后的非钻孔位置地层参数赋予有
限元模型。
11.本发明仅通过钻孔数据就能进行有限元地层建模，无需借助地质断面图，避免了随机场方法庞大的计算量，计算简单，使用方便。本发明在有限元计算中实现了多地层多参数建模，可对地层起伏工况进行计算，能更加准确的反应地层不均匀性对计算结果的影响。
12.步骤s1中，建立无地层参数的有限元模型的具体实现过程包括：
13.设置有限元模型的几何尺寸大小，完成有限元模型整体尺寸的几何建模，得到有限元整体几何模型；
14.对有限元整体几何模型进行网格划分，并对划分后的单元和节点进行编号；
15.导出所述单元和节点的对应关系以及各节点的坐标值。
16.本发明在几何建模时考虑地层起伏特性将各地层分层建出极易导致网格畸变问题，对模型整体尺寸进行几何建模，建模时不在几何建模阶段考虑地层起伏特性，是在整体模型上进行网格划分，通过对均匀网格赋予不同地层参数的手段既能体现地层起伏特性，又能解决几何划分地层所带来的网格畸变问题，极大地减小了计算结果误差，提高了计算精度。
17.优选地，为便于后续计算，所述单元和节点的对应关系存储于矩阵a中，所述矩阵a第一列为单元编号，其余列为节点编号；矩阵a每一行中第一个数字为单元编号，其余数字均为该单元所对应的节点编号。
18.优选地，为便于后续计算，所述各节点的坐标值存储于矩阵b中，矩阵b第一列为节点编号，其余列为坐标值；其中矩阵b每一行中第一个数字为节点编号，其余数字均为该节点编号对应的坐标值。
19.步骤s2中，将钻孔数据离散化，对所述有限元模型赋予相应参数的具体实现过程包括：
20.计算所述有限元模型中各单元中心点的坐标值；
21.根据所述各单元中心点的坐标值，首先查找平面坐标(除埋深坐标外其他坐标信息)等于各钻孔平面坐标的单元，若没有与某钻孔平面坐标相同的单元，则查找与该钻孔平面坐标距离最近的单元，将查找的所有单元组成矩阵e；
22.对比钻孔内的土层埋深和单元的埋深坐标信息，将钻孔相应埋深的参数赋予所述矩阵e。
23.经过步骤s2的处理后，有限元模型已从原来的无参数状态转变为带有钻孔参数的有限元模型，即模型中位于钻孔位置的单元已经被赋予了参数。
24.步骤s3的具体实现过程包括：
25.将各单元中心点的坐标值存储于矩阵d中；
26.将所述矩阵d与矩阵e做差，得到非钻孔位置单元的中心坐标集合f。
27.本发明中，步骤s4的具体实现过程包括：
28.将非钻孔位置第i个单元参数f
ij
与钻孔位置每个单元的参数相减，挑选出最小的差值，将该差值最小的钻孔位置单元参数作为非钻孔位置第i个单元参数，代替f
ij
；
29.遍历完所有非钻孔位置单元，完成全部非钻孔位置单元参数的更新；更新后的每个非钻孔位置单元的参数都是钻孔位置单元参数的一种，完成单元参数分类，即实现地层分类；
30.将更新好的单元参数按照单元编号赋予到有限元模型中，得到全地层参数的有限元模型；
31.其中，e
mj
为钻孔位置单元矩阵e中第m个单元的第j个力学参数；n为钻孔位置单元矩阵e的总个数；p
im
为计算非钻孔位置第i个单元参数时钻孔位置单元矩阵e中第m个单元参数所占权重，dx
im
、dy
im
分别为钻孔位置第m个单元和非钻孔位置第i个单元之间的水平向和竖向距离，dx
il
、dy
il
分别为钻孔位置第l个单元和非钻孔位置第i个单元之间的水平向和竖向距离，k为比例系数，n为钻孔位置单元矩阵e中的单元个数。
32.步骤s4在计算权重时对固有反距离加权平均方法进行了改进，引入了比例系数k，改进后的方法可以充分考虑地层力学参数在不同方向上的相关性。反距离加权平均法是地质中的概念，属于插值方法中的一种，即通过几个已知点的高程来计算未知点的高程。因高程信息无明显的方向相关性，所以反距离加权平均法效果尚可。但地层力学参数信息不同高程信息，大部分的地层都是沿水平方向分布的，具有明显的方向相关性。本发明中，通过改变数值k的大小，即可反应地层的不同方向的相关性(即地层主要沿水平分布还是竖直分布)。其中k等于1时表示x方向与y方向距离影响强度相同。k值大于1则代表单元间水平向距离对参数的作用强于竖向距离；k值小于1则代表单元间竖向距离对参数的影响大于水平向距离。由于一般情况下地层水平方向相关性大于竖直方向相关性(即地层在水平方向连续性强，大部分沿水平方向分布)，因此k值一般取大于1的数。
33.当有限元模型为三维模型时，p
im
通过下式计算：
[0034][0035]
通过上式，即可根据钻孔数据建立三维有限元力学参数模型，拓展了本发明方法的适用范围。
[0036]
作为一个发明构思，本发明还提供了一种终端设备，包括存储器、处理器及存储在存储器上的计算机程序；所述处理器执行所述计算机程序，以实现本发明上述方法的步骤。
[0037]
作为一个发明构思，本发明还提供了一种计算机可读存储介质，其上存储有计算机程序/指令；所述计算机程序/指令被处理器执行时实现本发明上述方法的步骤。
[0038]
与现有技术相比，本发明所具有的有益效果为：
[0039]
1)本发明在有限元计算中实现了多地层多参数建模，可对地层起伏工况进行计算，能更加准确的反应地层不均匀性对计算结果的影响；
[0040]
2)本发明仅通过钻孔数据就能进行有限元地层建模，无需借助地质断面图，使用更加方便；
[0041]
3)本发明的方法通过对均匀网格赋予不同地层参数的方案来体现地层起伏特性，解决了有限元模型中对地层直接进行几何划分所带来的网格畸变问题；
[0042]
4)本发明的方法所得地层参数为确定值，对于单个单元来说参数是唯一的，因此仅需计算一次就能得到稳定结果。相较于借助随机场方法需要大量计算得到统计学稳定值来说，计算成本低；
[0043]
5)本发明避免了直接在有限元模型进行地层几何划分所带来的网格畸变等问题，可根据已有钻孔数据在网格均一的有限元模型中生成不同地层参数，以此反映不同地层起伏变化规律，从而实现带倾斜地层的有限元数值模拟计算。
附图说明
[0044]
图1为本发明实施例1的方法步骤流程图；
[0045]
图2为本发明实施例1有限元模型网格划分及单元节点编号示意图。
[0046]
图3(a)和图3(b)为本发明实施例1带钻孔参数的有限元模型，其中图3(a)为有限元模型弹模参数的分布情况，图3(b)为有限元模型泊松比参数分布情况。
[0047]
图4(a)和图4(b)为本发明实施例1土体初始参数的有限元模型，其中图4(a)为有限元模型弹模参数的分布情况，图4(b)为有限元模型泊松比参数分布情况。
[0048]
图5(a)和图5(b)为本发明实施例1土体最终参数的有限元模型，其中图5(a)为有限元模型弹模参数的分布情况，图5(b)为有限元模型泊松比参数分布情况。
具体实施方式
[0049]
下面结合附图和实例对本发明的技术方案进行具体说明。需要说明的是，为了简单方便的对本发明进行阐述，选取了二维工程实例进行说明，并不说明本发明仅可用于二维有限元建模，采用同样的方法可以用与三维有限元建模中。
[0050]
实施例1
[0051]
针对某具体工程，需要建立长50m，埋深50m的二维有限元模型。在该范围内有两个钻孔数据，其中钻孔1的信息如下：钻孔位于长度方向9m处，埋深50m范围内岩土体共分为3层，埋深由上而下分别为20m，10m，20m，弹性模量分别为100mpa，20mpa，100mpa，泊松比分别为0.3,0.2,0.3。钻孔2的信息如下：钻孔位于长度方向43m处，埋深50m范围内岩土体共分为3层，埋深由上而下分别为10m，30m，20m，弹性模量分别为100mpa，20mpa，100mpa，泊松比分别为0.3,0.2,0.3。
[0052]
采用本发明实施例1的一种基于钻孔数据的地层力学参数有限元建模方法，按照以下步骤进行处理：
[0053]
步骤1)建立无地层参数的有限元模型：
[0054]
根据研究问题建立长50m(x方向)，埋深50m(y方向)的几何模型；随后对几何模型进行网格划分，并对单元(网格单元)和节点进行编号；为了简单的呈现本发明实施例1的方法，网格划分尺寸选择2.5m。网格划分及完成单元节点编号如图2所示，共划分了400个网格，得到了无地层参数的有限元模型。
[0055]
导出单元及节点信息，具体包括单元与节点的对应关系以及各节点的坐标值；
[0056]
为方便后续数据处理，将单元与节点对应关系存储于矩阵a中，其中矩阵a第一列为单元编号，其余列为节点编号。其中矩阵a每一行中第一数字为单元编号，其余数字均为该单元所对应的节点编号。本实施例中矩阵a如下：
[0057]
[0058]
优选地，为方便后续数据处理，各节点坐标值存储于矩阵b中，其中矩阵b第一列为节点编号，其余列为坐标值；其中矩阵b每一行中第一数字为节点编号，其余数字均为该节点编号对应的坐标值。
[0059][0060]
步骤s2的工程现场数据采集的内容主要包括：
[0061]
各钻孔位置及编号，各钻孔内各土层的埋深及力学参数；
[0062]
优选地，为方便程序后续批量作业，所述各钻孔内各土层的埋深及力学参数信息可采用矩阵c进行表示，其中矩阵c第一列为各岩土层埋深厚度，矩阵c其余列为对应的各岩土层的力学参数，如重度、弹性模量、泊松比、粘聚力、内摩擦角等。
[0063]
本实施例共有两个钻孔，则用分别用两个矩阵c1,c2进行信息存储，如下：
[0064][0065][0066]
步骤s3的具体实现过程包括：
[0067]
计算步骤s1中有限元模型中各单元中心点的坐标值，具体通过步骤s1所述的矩阵a与矩阵b计算求得，其中在矩阵a中提取各单元对应的节点编号，然后根据节点编号到矩阵b中查找每个节点的坐标值，分别对单元对应所有节点的x,y,z坐标值分别求平均(二维模型中对x，y坐标值求平均)，即为各单元中心点坐标值，将单元编号与坐标值信息存储到矩阵d中；
[0068]
根据上述所获得的各单元的中心点坐标位置信息(矩阵d)，首先查找平面坐标(除埋深坐标外其他坐标信息)等于各钻孔平面坐标的单元，若没有与某钻孔平面坐标相同的单元，则查找与该钻孔平面坐标距离最近的单元，将查找的单元形成矩阵e。随后通过对比钻孔内的土层埋深和单元的埋深坐标信息，将钻孔相应埋深的参数赋予上述钻孔位置单元矩阵e，以此完成步骤s1所述的网格化的有限元模型中赋予钻孔位置单元地层力学参数。
[0069]
经过步骤s3的处理后，有限元模型已从原来的无参数状态转变为带有钻孔参数的有限元模型，即模型中位于钻孔位置的单元已经被赋予了参数，如图3(a)和图3(b)所示。
[0070]
步骤s4的具体实现过程包括：
[0071]
将步骤s3中所有单元中心坐标(矩阵d)与钻孔位置单元中心坐标(矩阵e)作差集，得到非钻孔位置单元的中心坐标集合f。
[0072]
其中非钻孔位置中心坐标集合f中第i个单元的第j个力学参数f
ij
通过下式取得：
[0073][0074]
上式中，e
mj
为钻孔位置单元矩阵e中第m个单元的第j个力学参数；
[0075]
n为单元矩阵e的总个数。
[0076]
p
im
为关于所求非钻孔位置第i个单元参数时钻孔位置单元矩阵e中第m个单元参数
所占权重，具体计算方法如下式：
[0077][0078]
上式中：
[0079]
下标i表示非钻孔位置单元集合的第i个单元；
[0080]
下标m、l分别表示钻孔位置单元矩阵e的第m个、第l个单元；
[0081]
dx，dy分别为两个单元之间的水平向和竖向距离；
[0082]
k为比例系数，反映单元间水平向距离与竖向距离对土体参数影响强度的比值，其中k等于1时表示x方向与y方向距离影响强度相同。由于本实例所处区域地层在水平向的相关性远大于竖向相关性，所以取k＝10。
[0083]
优选地，当有限元模型为三维模型时，p
im
通过下式计算：
[0084][0085]
上式中dx，dy表示两个单元间水平向的距离，dz表示两个单元间的竖向距离。
[0086]
经过步骤s4处理后，有限元模型的非钻孔位置的单元就有了土体初步参数，该初步参数一般与钻孔原始参数有所差别，若将初始参数直接赋予有限元模型，则会产生图4(a)和图4(b)的结果。
[0087]
步骤s5的具体实现过程包括以下步骤：
[0088]
将非钻孔位置第i个单元参数f
ij
与钻孔位置每个单元的参数进行相减，挑选差值绝对值最小的，将该钻孔位置单元参数作为非钻孔位置第i个单元参数以代替f
ij
。采用该方法遍历完所有非钻孔位置单元，以此完成全部非钻孔位置单元参数的更改。更改后的每个非钻孔位置单元的参数都是钻孔位置单元参数的一种，以此完成单元参数分类，即实现了地层分类。
[0089]
将更新好的单元参数按照单元编号赋予到有限元模型中，以此实现全地层参数的有限元模型，全地层参数有限元模型如图5(a)和图5(b)所示。
[0090]
实施例2
[0091]
本发明实施例3提供一种对应上述实施例1的终端设备，终端设备可以是用于客户端的处理设备，例如手机、笔记本电脑、平板电脑、台式机电脑等，以执行上述实施例的方法。
[0092]
本实施例的终端设备包括存储器、处理器及存储在存储器上的计算机程序；处理器执行存储器上的计算机程序，以实现上述实施例1方法的步骤。
[0093]
在一些实现中，存储器可以是高速随机存取存储器(ram：random access memory)，也可能还包括非不稳定的存储器(non-volatile memory)，例如至少一个磁盘存储器。
[0094]
在另一些实现中，处理器可以为中央处理器(cpu)、数字信号处理器(dsp)等各种类型通用处理器，在此不做限定。
[0095]
实施例3
[0096]
本发明实施例3提供了一种对应上述实施例的计算机可读存储介质，其上存储有计算机程序/指令。计算机程序/指令被处理器执行时，实现上述实施例1方法的步骤。
[0097]
计算机可读存储介质可以是保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以是但不限于电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意组合。
[0098]
本领域内的技术人员应明白，本技术的实施例可提供为方法、系统、或计算机程序产品。因此，本技术可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本技术可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。本技术实施例中的方案可以采用各种计算机语言实现，例如，面向对象的程序设计语言java和直译式脚本语言javascript等。
[0099]
本技术是参照根据本技术实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0100]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0101]
尽管已描述了本技术的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本技术范围的所有变更和修改。
[0102]
显然，本领域的技术人员可以对本技术进行各种改动和变型而不脱离本技术的精神和范围。这样，倘若本技术的这些修改和变型属于本技术权利要求及其等同技术的范围之内，则本技术也意图包含这些改动和变型在内。