CN110826276B - 一种基于有限元的三维随机渗流场模拟方法及装置 - Google Patents

Abstract

本发明公开了一种基于有限元的三维随机渗流场模拟方法及装置，属于岩土工程随机渗流场模拟领域，基于有限元网格积分点坐标，将材料的渗透系数与温度场通过UTEMP子程序间接建立联系，能够模拟实现在三维空间上不均匀的三维随机渗流场；为满足材料渗透等级的要求，将渗透系数视为服从均匀统计分布的随机变量，同时可调整自相关函数中的相关长度来反映材料的渗透性质在空间各方向上的相关性。将三维随机渗流场的生成与有限元计算高效结合，可运用于任意复杂形状的三维随机渗流场模型，有效提高了计算能力和效率。有限元计算结果能够真实反映在材料渗透系数不均匀的情况下渗流状态的变化规律，有助于进一步明确模型总渗流量的不确定性来源。

Description

一种基于有限元的三维随机渗流场模拟方法及装置

技术领域

本发明属于岩土工程随机渗流场模拟技术领域，更具体地，涉及一种基于有限元的三维随机渗流场模拟方法及装置。

背景技术

水在孔隙介质中的流动称为渗流，水流所占有的空间区域称为渗流场。在岩土工程和水利水电工程中，因水的渗流问题导致的工程事故屡见不鲜。随着时间的推移，山体滑坡、地下隧道坍塌、深基坑垮塌、土石坝溃坝、城市地面下陷等工程事故近年来持续不断，工程的安全问题引起了科研人员和工程师的高度重视。究其根本原因，并非本身的强度不够，很大程度上是由于水流在岩土体结构内的渗流而引起的。由于地质构造的复杂性，岩土体在自然状态下呈现明显的非均质性和各向异性。岩土介质在不同区域的孔隙率存在较大的差异，对过水能力造成一定的影响，使得渗透系数在空间上存在较大的空间变异性。岩土介质渗透系数的不确定性进而导致土体内部的渗流场具有明显的随机性，为工程风险评估带来较大的不确定性。因此，如何真实反演随机渗流场的运动状态具有重要的工程实践和指导意义。

在现有的研究当中，随机场被公认为是表征材料空间不均匀性的有效途径。由于岩土材料的各项参数一般为非负数，为保证随机场产生的值为正值，通常假定渗透系数服从对数正态统计分布。在相近区域，岩土体材料的性质比较接近。根据规范《水利水电工程地质勘察规范》(GB50487-2008)，岩土体渗透性有明确的渗透分级标准，同一类岩土体的渗透系数一般在1～2个数量级范围内。若采用对数正态分布进行模拟，产生的渗透系数随机值将有一部分远远超过这个范围，从某种意义上来说，已经改变了土体的真实性质，这对渗流分析结果无疑将造成较大的影响。因此，在利用随机场来模拟渗透系数的空间变异性进行随机渗流分析时，有必要采取更合适的统计分布，确保渗透系数的随机值在同一类土体渗透系数变化的范围内，以获得更接近实际的数值模拟结果。

此外，在岩土工程实践中，工程中的渗流问题实际上属于复杂的三维空间问题。由于数值方法和技术所限，在一般的土石坝、基坑、隧洞等渗流场分析中通常将其简化为平面应变问题进行研究。基于二维随机渗流场得到的结果能在一定程度上反映渗流的运动状态，但与空间上的实际情况还存在一定的差距。渗流场各个区域的土质条件、边界条件在三维空间各个方向都是动态变化的，二维渗流模型不能充分体现这些特性，使得模拟结果较为粗略。因此，亟需发展一种能够模拟三维随机渗流场的方法，才可能达到对实际复杂情况的有效模拟。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提出了一种基于有限元的三维随机渗流场模拟方法及装置，由此解决现有渗流场模拟方法不能有效反映渗流场在三维空间上分布的不均匀性的技术问题。

为实现上述目的，按照本发明的一个方面，提供了一种基于有限元的三维随机渗流场模拟方法，包括：

(1)根据ABAQUS的CAE方式创建三维模型，将所述三维模型导出为INP文件，在所述INP文件中将材料的渗透系数与温度场变量进行关联，加入场变量语句，设置单元输出变量；

(2)利用ABAQUS自带的UTEMP子程序进行参数传递，确定能够进行参数传递的输入输出变量；

(3)假设渗透系数的值服从均匀统计分布，确定渗透系数变化范围的最小值和最大值，在UTEMP子程序自定义区域中利用Fortran语言编写程序，使用修正线性估计法生成三维随机场；

(4)调整所述三维随机场中自相关函数中渗透系数在x、y和z三个方向上的相关长度，以反映材料参数在各方向上的相关性；

(5)在不启动ABAQUS软件的情况下，通过DOS界面输入命令调用连接所述INP文件和所述UTEMP子程序，在后台进行有限元计算；

(6)打开.odb文件，输出所述三维模型的三维随机渗流场云图，以观察渗流速度在各区域的变化情况。

优选地，步骤(1)包括：

启动ABAQUS软件，根据CAE界面方式创建三维模型，设置流体渗透/应力耦合Soils分析；

采用C3D8P孔压单元类型将所述三维模型自由划分为六面体单元，导出含所述三维模型各积分点坐标的INP文件；

在所述INP文件中修改渗透系数的设置时，将渗透系数与温度场间接建立关系，在所述INP文件中相应地方添加与温度场变量相关的语句，并将TEMP设置为单元输出变量。

优选地，步骤(2)包括：

使用ABAQUS自带的与温度有关的UTEMP子程序中的传入变量COORDS和定义变量TEMP进行参数传递，在各所述积分点坐标处输出的温度值在数值上等价于渗透系数，能够间接实现渗透系数的随机性。

优选地，步骤(3)包括：

假设渗透系数的取值服从均匀统计分布，并确定所述三维模型中渗透系数变化范围的最小值和最大值；

在UTEMP子程序自定义区域中利用Fortran语言编写程序，使用修正线性估计法生成三维均匀分布随机场。

优选地，步骤(4)包括：

利用随机场表征渗透系数的空间不均匀性时，通过调整自相关函数中三个方向的相关长度，以反映材料的渗透系数在各方向上的相关性。

优选地，步骤(5)包括：

在不启动ABAQUS的情况下，打开DOS界面并将路径调至当前工作目录，输入命令连接所述INP文件和所述UTEMP子程序，在后台调用ABAQUS程序进行有限元计算，以提高计算效率和节约计算资源。

优选地，步骤(6)包括：

在工作目录下找到并打开.odb结果文件，在结果输出选项中，选择TEMP作为输出变量，显示所述三维模型的随机渗流场云图，并观察渗流速度在各区域的变化情况。

按照本发明的另一个方面，提供了一种基于有限元的三维随机渗流场模拟装置，包括：

模型构建模块，用于根据ABAQUS的CAE方式创建三维模型，将所述三维模型导出为INP文件，在所述INP文件中将材料的渗透系数与温度场变量进行关联，加入场变量语句，设置单元输出变量；

参数传递模块，用于利用ABAQUS自带的UTEMP子程序进行参数传递，确定能够进行参数传递的输入输出变量；

随机场生成模块，用于假设渗透系数的值服从均匀统计分布，确定渗透系数变化范围的最小值和最大值，在UTEMP子程序自定义区域中利用Fortran语言编写程序，使用修正线性估计法生成三维随机场，并调整所述三维随机场中自相关函数中渗透系数在x、y和z三个方向上的相关长度，以反映材料参数在各方向上的相关性；

计算模块，用于在不启动ABAQUS软件的情况下，通过DOS界面输入命令调用连接所述INP文件和所述UTEMP子程序，在后台进行有限元计算；

输出模块，用于打开.odb文件，输出所述三维模型的三维随机渗流场云图，以观察渗流速度在各区域的变化情况。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

(1)本发明将二维随机渗流场拓展至三维，有效克服了传统平面应变模型在边界条件设置方面的缺陷。三维随机渗流模型允许渗透系数在垂直平面的方向上发生变化，真正实现了渗透系数在三维空间上的不均匀性；

(2)本发明充分利用ABAQUS中温度子程序UTEMP的优势，将渗透系数与温度场在数值上建立等效关系，间接反映了渗透系数在三维空间的随机特性；

(3)由于同一区域渗透材料的渗透等级一般在1～2个数量级范围内，若采用对数正态分布，产生的随机场将有一部分值远远超过既定的数量级范围，不符合事实的要求。本发明视渗透系数服从均匀统计分布，打破了传统方法假定渗透系数服从对数正态分布的界限；

(4)本发明将三维随机渗流场与有限元高效结合，有限元结果能够真实反映渗流状态随渗透系数空间的不均匀性而发生的动态变化，与传统的确定性分析结果形成鲜明对比；

(5)本发明所提方法能够运用到任意复杂的三维随机渗流模型中，将三维随机场的高效生成与有限元的高度非线性计算融为一体，大大提高了复杂三维模型进行有限元计算的能力和效率，节约了计算时间和资源成本。

附图说明

图1是本发明实施例提供的一种基于有限元的三维随机渗流场模拟方法的流程图；

图2是本发明实施例提供的一种透水地基上的土石坝三维有限元模型的示意图，其中，上游孔压边界为65kPa，下游孔压边界为0kPa；

图3是本发明实施例提供的一种土石坝坝基三维随机渗流场的典型示意图，其中，坝体的渗透系数为确定值1×10^-5m/s，防渗墙为不透水材料；

图4是本发明实施例提供的一种图3中区域A的局部放大图，其中，(a)为实施例中土石坝三维渗流模型确定性分析在断面Z＝10m处的渗流状态，(b)为实施例中土石坝三维渗流模型确定性分析在断面Z＝20m处的渗流状态，(c)为实施例中土石坝三维渗流模型确定性分析在断面Z＝30m处的渗流状态，其中，箭头代表流速的大小和方向；

图5是本发明实施例提供的一种图3中区域A的局部放大图，其中，(a)为实施例中土石坝三维渗流模型随机分析在断面Z＝10m处的渗流状态的对比图，(b)为实施例中土石坝三维渗流模型随机分析在断面Z＝20m处的渗流状态的对比图，(c)为实施例中土石坝三维渗流模型随机分析在断面Z＝30m处的渗流状态的对比图，其中，箭头代表流速的大小和方向。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明通过一个无限深透水地基上的土石坝坝基三维随机渗流工程实施例，并结合附图，对本发明提出的技术方案做进一步阐述。

如图1所示是本发明实施例提供的一种基于有限元的三维随机渗流场模拟方法，包括以下步骤：

S1：启动ABAQUS软件，如图2所示，根据CAE方式拉伸创建三维有限元模型，采用流体渗透/应力耦合Soils分析步。在Mesh模块中，采用C3D8P孔压单元类型将三维模型自由划分为六面体单元，导出INP文件，内含模型各积分点的坐标。在INP文件中修改渗透系数的设置时，由于在CAE中没有能够与渗透系数直接相关联的场变量，本发明实施例通过将渗透系数与温度场间接建立等价关系。在INP中相应的地方添加温度场变量语句，并将TEMP设置为单元输出变量。插入语句的具体内容如下。

1)将坝基的渗透系数与温度进行关联(第一列与第三列数据需保持一致)：

*Permeability,specific＝10.

5e-09,0,5e-09,

0.005,0,0.005,

0.05,0,0.05,

5,0,5,

2)在预定应场PREDEFINED FIELDS中设置初始温度为0，其中Base为模型中坝基的单元集合：

*Initial Conditions,TYPE＝TEMPERATURE

Base,0

3)在分析步中边界条件之后加入调用UTEMP子程序语句：

*TEMPERATURE,USER

Base,

4)设置单元温度场变量输出：

*Element Output,directions＝YES

TEMP

S2：在本发明实施例中，间接使用ABAQUS自带的温度子程序UTEMP进行参数传递，主要通过传入变量COORDS和定义变量TEMP传递参数，从而实现坝基渗透系数的随机性。在各积分点坐标处输出的温度值在数值上等价于渗透系数，然后通过插值得到模型渗透系数的云图。UTEMP子程序中固定的格式如下：

SUBROUTINE UTEMP(TEMP,NSECPT,KSTEP,KINC,TIME,NODE,COORDS)

INCLUDE'ABA_PARAM.INC'

DIMENSION TEMP(NSECPT),TIME(2),COORDS(3)

用户代码定义U

RETURN

END

由于同一区域渗透材料的渗透等级一般在1～2个数量级范围内，为基本满足这一事实要求，在本发明实施例中，假设渗透系数的取值服从均匀统计分布，并确定模型中渗透系数变化范围的最小值和最大值。在UTEMP子程序自定义区域中利用Fortran语言编写程序，使用修正线性估计法生成三维均匀分布随机场。

S3：由于同一区域渗透材料的渗透等级一般在1～2个数量级范围内，为基本满足这一事实要求，在本发明实施例中，假设渗透系数的取值服从均匀统计分布，并确定模型中渗透系数变化范围的最小值和最大值。对于本发明实施例，坝基地层主要为粉细砂透水层，渗透性等级为弱透水，渗透系数K的变化范围约为10^-6到10^-4。本发明实施例将此范围的极值分别作为服从均匀统计分布的渗透系数取值范围的最小值k_min和最大值k_max。在子程序UTEMP用户代码定义区域中利用Fortran语言编写程序，使用修正线性估计法生成三维均匀分布随机场，在积分点处生成的随机值在数值上等于材料的渗透系数。

1)修正线性估计法的公式为：

式(1)中，x为总体坐标系下的积分点坐标(x,y,z)；

N_i(x)为各点处的形函数，其中s、t、z′分别为自然坐标系下的相对位置坐标：

N₁＝[(1+s)×(1+t)×(1+z′)]/8

N₂＝[(1-s)×(1+t)×(1+z′)]/8

N₃＝[(1-s)×(1-t)×(1+z′)]/8

N₄＝[(1+s)×(1-t)×(1+z′)]/8

N₅＝[(1+s)×(1+t)×(1-z′)]/8

N₆＝[(1-s)×(1+t)×(1-z′)]/8

N₇＝[(1-s)×(1-t)×(1-z′)]/8

N₈＝[(1+s)×(1-t)×(1-z′)]/8；

f_i为自然坐标系下八角点处的随机值。

2)均匀分布随机场的公式为：

temp(1)＝k_min+(k_max-k_min)×f(x) (2)

式(2)中，temp(1)为积分点处输出的温度值，在数值上等于渗透系数；

k_min、k_max分别为渗透系数取值范围的最小值和最大值。

S4：自然岩土体的渗透系数在空间一定范围内往往具有较强的相关性，利用随机场表征渗透系数的空间不均匀性时，通过调整自相关函数中三个方向的相关长度可以有效反映材料参数在各方向上的相关性。岩土工程中，较常用的相关函数为平方指数型自相关函数，表达形式如下：

式(3)中：ρ(x,y,z)为自相关函数；

Δx、Δy和Δz分别为空间两点之间在x、y和z方向的相对距离；

θ_x、θ_y和θ_z分别为渗透系数在x、y和z方向的相关长度，在本实施例中，三者分别为20m、2m和20m。

S5：为提高计算效率，在不启动ABAQUS软件的情况下，将INP文件和子程序放在同一工作目录下，打开DOS界面并将路径调至当前工作目录，输入命令连接INP文件和UTEMP子程序，在后台调用ABAQUS程序进行三维有限元计算。具体调用的格式为：

abaqus job＝“job-name”user＝“sub-name”int

其中：job-name为INP文件名；sub-name为子程序名，不需要加后缀。

S6：在工作目录下找到并打开.odb结果文件，在结果输出选项中，选择TEMP作为输出变量，显示土石坝坝基三维随机渗流场云图，如图3所示，并观察渗流速度在各区域的变化情况，如图5所示，其中，图5中(a)为实施例中土石坝三维渗流模型随机分析在断面Z＝10m处的渗流状态的对比图，图5中(b)为实施例中土石坝三维渗流模型随机分析在断面Z＝20m处的渗流状态的对比图，图5中(c)为实施例中土石坝三维渗流模型随机分析在断面Z＝30m处的渗流状态的对比图，可见流速在Z方向上各断面的状况具有明显差异，体现了渗流在三维空间上的各向异性，反映了真实的渗流情况。与土石坝三维渗流模型的确定性分析结果形成鲜明对比，如图4所示，其中，图4中(a)为实施例中土石坝三维渗流模型确定性分析在断面Z＝10m处的渗流状态，图4中(b)为实施例中土石坝三维渗流模型确定性分析在断面Z＝20m处的渗流状态，图4中(c)为实施例中土石坝三维渗流模型确定性分析在断面Z＝30m处的渗流状态，其中，箭头代表流速的大小和方向，进一步说明了本发明所提方法的有效性和实用性。

在本发明的另一实施例中，还提供了一种基于有限元的三维随机渗流场模拟装置，包括：

模型构建模块，用于根据ABAQUS的CAE方式创建三维模型，将所述三维模型导出为INP文件，在所述INP文件中将材料的渗透系数与温度场变量进行关联，加入场变量语句，设置单元输出变量；

参数传递模块，用于利用ABAQUS自带的UTEMP子程序进行参数传递，确定能够进行参数传递的输入输出变量；

随机场生成模块，用于假设渗透系数的值服从均匀统计分布，确定所述渗透系数变化范围的最小值和最大值，在UTEMP子程序自定义区域中利用Fortran语言编写程序，使用修正线性估计法生成三维随机场，并调整所述三维随机场中自相关函数中渗透系数在x、y和z三个方向上的相关长度，以反映材料参数在各方向上的相关性；

计算模块，用于在不启动ABAQUS软件的情况下，通过DOS界面输入命令调用连接所述INP文件和所述UTEMP子程序，在后台进行有限元计算；

输出模块，用于打开.odb文件，输出模型的三维随机渗流场云图，观察渗流速度在各区域的变化情况。

其中，各模块的具体实施方式可以参考方法实施例的描述，本发明实施例将不再复述。

本发明充分利用大型通用有限元软件ABAQUS的平台优势，基于有限元网格积分点坐标，将材料的渗透系数与温度场通过UTEMP子程序间接建立联系，能够模拟实现在三维空间上不均匀的三维随机渗流场，有效解决了传统平面应变模型中边界条件对渗流的影响作用，突破了二维随机渗流场模型和三维确定性模型的局限性。为满足材料渗透等级的要求，本发明通过将渗透系数视为服从均匀统计分布的随机变量，打破了传统方法中对数正态分布的界限，同时可调整自相关函数中的相关长度来反映材料的渗透性质在空间各方向上的相关性。本发明将三维随机渗流场的生成与有限元计算高效结合，可运用于任意复杂形状的三维随机渗流场模型，有效提高了计算能力和效率，节约了计算时间和资源。有限元计算结果能够真实反映在材料渗透系数不均匀的情况下渗流状态的变化规律，有助于进一步明确模型总渗流量的不确定性来源，为对大型三维复杂模型的随机渗流分析进行有限元计算提供了新的设计途径。

需要指出，根据实施的需要，可将本申请中描述的各个步骤/部件拆分为更多步骤/部件，也可将两个或多个步骤/部件或者步骤/部件的部分操作组合成新的步骤/部件，以实现本发明的目的。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种基于有限元的三维随机渗流场模拟方法，其特征在于，包括：

(1)根据ABAQUS的CAE方式创建三维模型，将所述三维模型导出为INP文件，在所述INP文件中将材料的渗透系数与温度场变量进行关联，加入场变量语句，设置单元输出变量；

(2)利用ABAQUS自带的UTEMP子程序进行参数传递，确定能够进行参数传递的输入输出变量；

(3)假设渗透系数的值服从均匀统计分布，确定渗透系数变化范围的最小值和最大值，在UTEMP子程序自定义区域中利用Fortran语言编写程序，使用修正线性估计法生成三维随机场；

(4)调整所述三维随机场中自相关函数中渗透系数在x、y和z三个方向上的相关长度，以反映材料参数在各方向上的相关性；

(5)在不启动ABAQUS软件的情况下，通过DOS界面输入命令调用连接所述INP文件和所述UTEMP子程序，在后台进行有限元计算；

(6)打开.odb文件，输出所述三维模型的三维随机渗流场云图，以观察渗流速度在各区域的变化情况。

2.根据权利要求1所述的方法，其特征在于，步骤(1)包括：

启动ABAQUS软件，根据CAE界面方式创建三维模型，设置流体渗透/应力耦合Soils分析；

采用C3D8P孔压单元类型将所述三维模型自由划分为六面体单元，导出含所述三维模型各积分点坐标的INP文件；

在所述INP文件中修改渗透系数的设置时，将渗透系数与温度场间接建立关系，在所述INP文件中相应地方添加与温度场变量相关的语句，并将TEMP设置为单元输出变量。

3.根据权利要求2所述的方法，其特征在于，步骤(2)包括：

使用ABAQUS自带的与温度有关的UTEMP子程序中的传入变量COORDS和定义变量TEMP进行参数传递，在各所述积分点坐标处输出的温度值在数值上等价于渗透系数，能够间接实现渗透系数的随机性。

4.根据权利要求3所述的方法，其特征在于，步骤(3)包括：

假设渗透系数的取值服从均匀统计分布，并确定所述三维模型中渗透系数变化范围的最小值和最大值；

在UTEMP子程序自定义区域中利用Fortran语言编写程序，使用修正线性估计法生成三维均匀分布随机场。

5.根据权利要求4所述的方法，其特征在于，步骤(4)包括：

利用随机场表征渗透系数的空间不均匀性时，通过调整自相关函数中三个方向的相关长度，以反映材料的渗透系数在各方向上的相关性。

6.根据权利要求1所述的方法，其特征在于，步骤(5)包括：

在不启动ABAQUS的情况下，打开DOS界面并将路径调至当前工作目录，输入命令连接所述INP文件和所述UTEMP子程序，在后台调用ABAQUS程序进行有限元计算，以提高计算效率和节约计算资源。

7.根据权利要求6所述的方法，其特征在于，步骤(6)包括：

在工作目录下找到并打开.odb结果文件，在结果输出选项中，选择TEMP作为输出变量，显示所述三维模型的随机渗流场云图，并观察渗流速度在各区域的变化情况。

8.一种基于有限元的三维随机渗流场模拟装置，其特征在于，包括：

模型构建模块，用于根据ABAQUS的CAE方式创建三维模型，将所述三维模型导出为INP文件，在所述INP文件中将材料的渗透系数与温度场变量进行关联，加入场变量语句，设置单元输出变量；

参数传递模块，用于利用ABAQUS自带的UTEMP子程序进行参数传递，确定能够进行参数传递的输入输出变量；

随机场生成模块，用于假设渗透系数的值服从均匀统计分布，确定渗透系数变化范围的最小值和最大值，在UTEMP子程序自定义区域中利用Fortran语言编写程序，使用修正线性估计法生成三维随机场，并调整所述三维随机场中自相关函数中渗透系数在x、y和z三个方向上的相关长度，以反映材料参数在各方向上的相关性；

计算模块，用于在不启动ABAQUS软件的情况下，通过DOS界面输入命令调用连接所述INP文件和所述UTEMP子程序，在后台进行有限元计算；

输出模块，用于打开.odb文件，输出所述三维模型的三维随机渗流场云图，以观察渗流速度在各区域的变化情况。

Images