CN105138738A - 一种三维渗透张量计算方法 - Google Patents

Abstract

本发明公开了一种三维渗透张量计算方法，在二维裂隙网络渗流分析的基础上，沿一定轴旋转二维剖面，计算不同旋转角度下，二维剖面各个方向的渗透系数，并绘制椭圆球体，最后用回归分析近似求得三维裂隙的渗透系数张量。由于在二维裂隙网络渗流计算中，采用一维的线单元模拟裂隙，有效地避免了三维裂隙网络模型中复杂的裂隙处理工作。在裂隙计算时，因为求解方程的维度由三维降低到二维，所以可以大幅度地提高有效裂隙的计算数量和效率。

Description

一种三维渗透张量计算方法

技术领域

本发明涉及一种三维渗透张量计算方法。

背景技术

经典渗流理论是以连续介质假定为基础的，而实际情况中，大多数介质都是离散的。为了能将离散介质适用于经典渗流理论，引入了表征单元体(REV)的概念，从而把实际上离散的多孔介质等效成由表征单元体组成的连续介质。当表征单元体尺寸与研究区域尺寸相比很小时，就可以用连续介质渗流理论进行分析。对于土体等孔隙介质，其表征单元体非常小，可以直接运用连续介质渗流理论。对于岩体，由于岩体中通常会有各种节理裂隙，在流体流经岩体时，主要是由裂隙通过的。裂隙和孔隙主要有两种差别：1.孔隙在三维方向上延伸，且各个方向上的差别不大，而裂隙主要在二维方向上延伸，其两个方向上的尺寸比第三个方向上的尺寸大很多。2.孔隙的延伸尺寸比研究域小得多，而裂隙的延伸尺度可以达到整个研究域。

一直以来，工程师在计算中把岩体看成和土体一样的多孔介质，应用经典渗流理论进行分析。而1959年法国Malpasset拱坝和1967年Vajiont拱坝事故的发生，使人们认识到岩体渗流不能和土体一样直接采用连续介质渗流理论分析。

如今，随着裂隙性油气田的开发、高放射性核废料的深埋处理和对拱坝岩体的开发，使得裂隙岩体渗流的研究逐渐成为了热点问题。目前对裂隙岩体进行渗流计算时，主要有：等效连续介质模型、离散裂隙网络模型和双重介质模型。但是由于后两者都有难以解决的问题而影响了模型的应用：在双重介质模型中，裂隙网络和连续介质岩块耦合模型中水交换量公式的准确性直接影响着模型的计算准确性，但是该模型中水交换量难以准确确定；离散裂隙网络模型认为岩块不能导水而忽略岩块的导水性，但实际情况中，岩块的非饱和渗透系数和裂隙渗透系数相比，一般不能忽略，故当考虑裂隙岩体非饱和渗流时，选用离散裂隙模型并不合适；相比之下，等效连续介质模型分析裂隙岩体渗流，在理论和实际中均有非常多的经验和成熟和计算程序可供借鉴。同时，该模型只需知道裂隙的几何信息的统计值，而对其准确位置并不敏感，因此，在实际工程问题中多使用等效连续介质模型。渗透张量是等效连续介质渗流计算模型中的一个关键参数，它直接反应了裂隙岩体非均质各向异性的特点。近年来，许多学者和专家在等效连续介质渗流计算中，就如何准确确定渗透张量这一问题进行了大量研究分析。

在三维情况下，用等效连续介质模型分析裂隙岩体的渗流时，需要使用一个3阶的渗透张量来表达裂隙渗流的各向异性特性。许多学者采用裂隙网络渗流的方法来计算三维渗透系数张量。LongJCS(1985)提出圆盘裂隙网络模型。Dershowitz(1987)提出多边形裂隙网络模型。CacasMC(1990)提出圆形管道模型。这些方法的实质是将裂隙及其交线看成一个点，将裂隙中的渗流简化成水在这些点之间的直线运动。但是，采用这些方法进行渗流分析时，由于三维情况下裂隙网络的复杂性，需要花费大量精力在裂隙生成、筛选以及裂隙与裂隙或边界的联系这些方面上。当岩体内部裂隙密度高时，可能在20米的尺度范围内就能包含上万条对渗流有作用的裂隙，这时，对裂隙的生成、筛选以及相关计算工作是非常巨大的。

发明内容

为了解决目前在三维情况下分析裂隙岩体的渗流计算工作量过大的技术问题，本发明提供一种可有效提高有效裂隙的计算效率的三维渗透张量计算方法。

为了实现上述技术目的，本发明的技术方案是：

一种三维渗透张量计算方法，包括以下步骤：

步骤一：输入由统计资料确定的各组裂隙面服从密度分布参数；

步骤二：在三维坐标中假定一个二维裂隙生成域，将优势裂隙面产状进行视倾角转换，获得裂隙在该二维生成域中迹长和方向角所服从的密度分布参数；

步骤三：在二维生成域中根据步骤二得到的迹长和方向角所服从的密度分布参数，采用ment-carlo法生成裂隙，并在生成域中心位置假定分析域尺寸，分析域的尺寸不大于生成域的尺寸；

步骤四：在分析域的范围内，对裂隙进行筛选，通过裂隙网络渗流模型计算分析域二维面上与渗透水流方向垂直的边界节点的流量Q和水头H，计算分析域二维平面上水流方向上的渗透系数K：

$K=\frac{V}{J}=\frac{Q/L}{(H_1-H_2)/L}=\frac{Q}{\mathrm{\Δ}H}---\left(1\right)$

其中Q和△H为与渗透水流方向垂直的边界节点的流量和水头差；V为流速，J为水力梯度，L为分析域尺寸，H₁，H₂为分析域区域左右两边的水头；

步骤五：在步骤三的二维生成域中，将分析域旋转角度β后重复步骤四，直到旋转一周为止，获得各次旋转后的渗透系数k_g(β)；

步骤六：将各个方向上的渗透系数经过处理后拟合成椭圆，计算椭圆拟合质量RMS_norm：

${\mathrm{RMS}}_{norm}=\frac{2}{(K_1+K_2)}\sqrt{\frac{1}{N}\underset{1}{\overset{N}{\Σ}}{\[k_g\left(\mathrm{\β}\right)-k_n\left(\mathrm{\β}\right)\]}^2}---\left(2\right)$

式中：N代表步骤五中计算得到的渗透系数的个数，k_g(β)是由旋转角度β后通过裂隙网络模型计算得到的渗透系数；k_n(β)是拟合椭圆上相应方向角β上的渗透系数；

步骤七：如果RMS_norm≤0.2则认为椭圆质量好，其表征单元体REV为分析域尺寸；反之，增大分析域尺寸，直到满足条件，但分析域尺寸大于生成域尺寸时，认为REV不存在；

步骤八：以任一穿过生成域的直线为定轴旋转二维裂隙生成域，直至旋转360°，重复步骤二～七，新旋转角度后的二维生成域的初始分析域尺寸取各REV中的最大值；

步骤九：根据前面计算出的各个方向的渗透系数及其方向，组成下式:

Ax＝b(3)

其中：

$A=\left[\begin{array}{cccccc}{n}_{11}^2& n_{12}^2& n_{13}^2& 2n_{11}{n}_{12}& 2n_{12}{n}_{13}& 2n_{11}{n}_{13}\\ n_{21}^2& n_{22}^2& n_{23}^2& 2n_{21}{n}_{22}& 2n_{22}{n}_{23}& 2n_{21}{n}_{23}\\ ...& ...& ...& ...& ...& ...\\ n_{n1}^2& n_{n2}^2& n_{n3}^2& 2n_{n1}{n}_{n2}& 2n_{n2}{n}_{n3}& 2n_{n1}{n}_{n3}\end{array}\right]$

x＝(K₁₁,K₂₂,K₃₃,K₁₂,K₂₃,K₁₃)^T

b＝(K_g(n₁),K_g(n₂),...,K_g(n_n))^T

式中：n_n1,n_n2,n_n3为的渗透系数方向的单位矢量的在x,y,z轴上分量；

K₁₁,K₂₂,K₃₃,K₁₂,K₂₃,K₁₃为三维渗透张量的六个独立分量；

K_g(n_i)为第i个方向的渗透系数；

根据最小二乘法求解出渗透系数张量的六个分量x，获得三维渗透张量K：

$K=\left[\begin{array}{ccc}{K}_{11}& K_{12}& K_{13}\\ K_{12}& K_{22}& K_{23}\\ K_{13}& K_{23}& K_{33}\end{array}\right].$

所述的一种三维渗透张量计算方法，所述的步骤二中，迹长和方向角所服从的密度分布参数包括迹长的均值和标准差，方向角的均值和标准差。

所述的一种三维渗透张量计算方法，所述的步骤二中，二位裂隙生成域的尺寸为裂隙全局域尺寸的一半，裂隙全局域尺寸为裂隙迹长的均值的2-4倍。

所述的一种三维渗透张量计算方法，所述的步骤二中，所述的视倾角γ转换通过tanγ＝tanαcosω计算，其中γ为视倾角，α为真倾角，ω为视倾向即剖面方向与倾向之夹角。

本发明的技术效果在于，在二维裂隙网络渗流分析的基础上，沿一定轴旋转二维剖面，计算不同旋转角度下，二维剖面各个方向的渗透系数，并绘制椭圆球体，最后用回归分析近似求得三维的渗透系数张量裂隙。由于在二维裂隙网络渗流计算中，采用一维的线单元模拟裂隙，有效地避免了三维裂隙网络模型中复杂的裂隙处理工作。在裂隙计算时，因为求解方程的维度由三维降低到二维，所以可以大幅度地提高有效裂隙的计算数量和效率。

具体实施方式

本方法包括以下步骤：

步骤一：输入由统计资料确定的各组裂隙面服从密度分布参数；

步骤二：在三维坐标中假定一个二维裂隙生成域，该二位裂隙生成域的尺寸为裂隙全局域尺寸的一半，裂隙全局域尺寸为裂隙迹长的均值的2-4倍，将优势裂隙面产状进行视倾角转换，其中，真倾向与视倾向之间夹角在倾斜面上斜交走向线所引的任一直线均为视倾斜线(如在任一斜交岩层走向的露头断面上所见)，视倾斜线与其在水平参考面上的投影线的夹角，叫视倾角。这样可获得裂隙在该二维生成域中迹长和方向角所服从的密度分布参数，当密度函数服从正态分布时，相关参数包括迹长的均值和标准差，方向角的均值和标准差；标准差和均值为统计分析学中的公知含义；

步骤三：在二维生成域中根据步骤二得到的迹长和方向角所服从的密度分布参数，采用ment-carlo法生成裂隙，并在生成域中心位置假定分析域尺寸，分析域的尺寸不大于生成域的尺寸；

步骤四：在分析域的范围内，对裂隙进行筛选，由于分析域是一个矩形，有两边与水流方向相同，两边与水流方向垂直，这里计算与水流方向垂直的边界上的节点，通过裂隙网络渗流模型计算分析域二维面上与渗透水流方向垂直的边界节点的流量Q和水头H，计算分析域二维平面上水流方向上的渗透系数K：

$K=\frac{V}{J}=\frac{Q/L}{(H_1-H_2)/L}=\frac{Q}{\mathrm{\Δ}H}---\left(1\right)$

其中Q和△H为与渗透水流方向垂直的边界节点的流量和水头差；V为流速，J为水力梯度，L为分析域尺寸，H₁，H₂为分析域区域左右两边的水头；

步骤五：在步骤三的二维生成域中，将分析域旋转角度β后重复步骤四，直到旋转一周为止，获得各次旋转后的渗透系数k_g(β)；

步骤六：将各个方向上的渗透系数经过处理后拟合成椭圆，计算椭圆拟合质量RMS_norm：

${\mathrm{RMS}}_{norm}=\frac{2}{(K_1+K_2)}\sqrt{\frac{1}{N}\underset{1}{\overset{N}{\Σ}}{\[k_g\left(\mathrm{\β}\right)-k_n\left(\mathrm{\β}\right)\]}^2}---\left(2\right)$

式中：N代表步骤五中计算得到的渗透系数的个数，k_g(β)是由旋转角度β后通过裂隙网络模型计算得到的渗透系数；k_n(β)是拟合椭圆上相应方向角β上的渗透系数；

步骤七：如果RMS_norm≤0.2则认为椭圆质量好，其表征单元体REV为分析域尺寸；反之，增大分析域尺寸，直到满足条件，但分析域尺寸大于生成域尺寸时，认为REV不存在；

步骤八：以任一穿过生成域的直线为定轴旋转二维裂隙生成域，直至旋转360°，重复步骤二～七，新旋转角度后的二维生成域的初始分析域尺寸取各REV中的最大值；

步骤九：根据前面计算出的各个方向的渗透系数及其方向，组成下式:

Ax＝b(3)

其中：

$A=\left[\begin{array}{cccccc}{n}_{11}^2& n_{12}^2& n_{13}^2& 2n_{11}{n}_{12}& 2n_{12}{n}_{13}& 2n_{11}{n}_{13}\\ n_{21}^2& n_{22}^2& n_{23}^2& 2n_{21}{n}_{22}& 2n_{22}{n}_{23}& 2n_{21}{n}_{23}\\ ...& ...& ...& ...& ...& ...\\ n_{n1}^2& n_{n2}^2& n_{n3}^2& 2n_{n1}{n}_{n2}& 2n_{n2}{n}_{n3}& 2n_{n1}{n}_{n3}\end{array}\right]$

x＝(K₁₁,K₂₂,K₃₃,K₁₂,K₂₃,K₁₃)^T

b＝(K_g(n₁),K_g(n₂),...,K_g(n_n))^T

式中：n_n1,n_n2,n_n3为的渗透系数方向的单位矢量的在x,y,z轴上分量；

K₁₁,K₂₂,K₃₃,K₁₂,K₂₃,K₁₃为三维渗透张量的六个独立分量；

K_g(n_i)为第i个方向的渗透系数；

根据最小二乘法求解出渗透系数张量的六个分量x，获得三维渗透张量K：

$K=\left[\begin{array}{ccc}{K}_{11}& K_{12}& K_{13}\\ K_{12}& K_{22}& K_{23}\\ K_{13}& K_{23}& K_{33}\end{array}\right].$

下面给出具体实施例：

某边坡主要发育有4组裂隙，具体产状如下：

第1组：层理面，集中产状：N45°W，NE∠11°，西侧基本贯通性连通，延伸长度大于20m，发育间距0.5m左右，宽度0.2mm～1.0mm；

第2组集中产状：N20°E，NW∠82°，除层间结构面以外，本组结构面最为发育，发育间距0.4m～1.0m，发育长度15m～30m，表现为张扭性特征，多张开0.5mm～1.5mm，少数与层间结构面交汇处发育溶槽；

第3组集中产状：N41°W，NE∠53°，主要发育在东侧边坡，发育间距0.4m～1.0m，发育长度5m～10m，多张开0.4mm～1.2mm；

第4组集中产状：N48°E，NW∠78°，主要发育在东侧边坡，发育间距0.5m～0.8m，发育长度15m～25m，宽度0.4mm～1.0mm。

各裂隙产状统计规律如下表：

表1裂隙产状参数表

在三维模型中以正北方为正Y轴，正东方为正X轴，向上方向为正Z轴，初始裂隙生成域设为XOZ平面。根据对该边坡节理裂隙产状描述，计算出节理裂隙在旋转后生成域中的二维裂隙方向角如下表2。由于二维裂隙是三维裂隙在该二维平面上的交线，所以已知三维裂隙的产状统计规律后，可以认为二维裂隙产状服从同样的规律。

表2不同旋转角度下裂隙方向角

最终计算出来的渗透张量(m/s)为：

$K=\left[\begin{array}{ccc}9.08\×{10}^{-6}& -3.57\×{10}^{-6}& -1.01\×{10}^{-6}\\ -3.57\×{10}^{-6}& 7.16\×{10}^{-6}& 1.56\×{10}^{-6}\\ -1.01\×{10}^{-6}& 1.56\×{10}^{-6}& 1.34\×{10}^{-6}\end{array}\right]$

对渗透张量求特征值，得出对应的渗透主值(m/s)为：

$K=\left[\begin{array}{ccc}4.38\×{10}^{-6}& & \\ & 1.07\×{10}^{-5}& \\ & & 1.45\×{10}^{-5}\end{array}\right]$

渗透主值对应的方向向量为：

$\mathrm{\ψ}=\left[\begin{array}{ccc}0.5932& -0.7013& -0.3952\\ 0.8018& 0.4707& 0.3682\\ -0.0722& -0.5353& 0.8416\end{array}\right]$

为了验证该方法的有效性，本文拟采用裂隙测量法对其进行验算。由对裂隙产状描述的内容发现，大多数裂隙节理延伸长度在15m以上，可认为在REV尺寸内是贯通的，采用裂隙测量法计算得到的渗透主值(m/s)为：

$K^{\′}=\left[\begin{array}{ccc}5.98\×{10}^{-6}& & \\ & 1.12\×{10}^{-5}& \\ & & 1.54\×{10}^{-5}\end{array}\right]$

渗透主值对应的方向向量为：

${\mathrm{\ψ}}^{\′}=\left[\begin{array}{ccc}0.6402& -0.6139& -0.4618\\ 0.7661& 0.4662& 0.4426\\ -0.0563& -0.6369& 0.7689\end{array}\right]$

对比发现，上述两种方法计算结果比较接近，在一定程度上证明了该方法的合理性。但是用裂隙测量法计算的结果偏大，这是由于假定裂隙是贯通的，而非连续裂隙网络计算中生成的裂隙并不全是贯通的，从而导致裂隙测量法结果偏大。

Claims (4)

1.一种三维渗透张量计算方法，其特征在于，包括以下步骤：

步骤一：输入由统计资料确定的各组裂隙面服从密度分布参数；

步骤二：在三维坐标中假定一个二维裂隙生成域，将优势裂隙面产状进行视倾角转换，获得裂隙在该二维生成域中迹长和方向角所服从的密度分布参数；

步骤三：在二维生成域中根据步骤二得到的迹长和方向角所服从的密度分布参数，采用ment-carlo法生成裂隙，并在生成域中心位置假定分析域尺寸，分析域的尺寸不大于生成域的尺寸；

步骤四：在分析域的范围内，对裂隙进行筛选，通过裂隙网络渗流模型计算分析域二维面上与渗透水流方向垂直的边界节点的流量Q和水头H，计算分析域二维平面上水流方向上的渗透系数K：

$K=\frac{V}{J}=\frac{Q/L}{(H_1-H_2)/L}=\frac{Q}{\mathrm{\Δ}H}---\left(1\right)$

其中Q和ΔH为与渗透水流方向垂直的边界节点的流量和水头差；V为流速，J为水力梯度，L为分析域尺寸，H₁，H₂为分析域区域左右两边的水头；

步骤五：在步骤三的二维生成域中，将分析域旋转角度β后重复步骤四，直到旋转一周为止，获得各次旋转后的渗透系数k_g(β)；

步骤六：将各个方向上的渗透系数经过处理后拟合成椭圆，计算椭圆拟合质量RMS_norm：

${\mathrm{RMS}}_{norm}=\frac{2}{(K_1+K_2)}\sqrt{\frac{1}{N}\underset{1}{\overset{N}{\Σ}}{\[k_g\left(\mathrm{\β}\right)-k_n\left(\mathrm{\β}\right)\]}^2}---\left(2\right)$

式中：N代表步骤五中计算得到的渗透系数的个数，k_g(β)是由旋转角度β后通过裂隙网络模型计算得到的渗透系数；k_n(β)是拟合椭圆上相应方向角β上的渗透系数；

步骤七：如果RMS_norm≤0.2则认为椭圆质量好，其表征单元体REV为分析域尺寸；反之，增大分析域尺寸，直到满足条件，但分析域尺寸大于生成域尺寸时，认为REV不存在；

步骤八：以任一穿过生成域的直线为定轴旋转二维裂隙生成域，直至旋转360°，重复步骤二～七，新旋转角度后的二维生成域的初始分析域尺寸取各REV中的最大值；

步骤九：根据前面计算出的各个方向的渗透系数及其方向，组成下式:

Ax＝b(3)

其中：

$A=\left[\begin{array}{cccccc}{n}_{11}^2& n_{12}^2& n_{13}^2& 2n_{11}{n}_{12}& 2n_{12}{n}_{13}& 2n_{11}{n}_{13}\\ n_{21}^2& n_{22}^2& n_{23}^2& 2n_{21}{n}_{22}& 2n_{22}{n}_{23}& 2n_{21}{n}_{23}\\ ...& ...& ...& ...& ...& ...\\ n_{n1}^2& n_{n2}^2& n_{n3}^2& 2n_{n1}{n}_{n2}& 2n_{n2}{n}_{n3}& 2n_{n1}{n}_{n3}\end{array}\right]$

x＝(K₁₁,K₂₂,K₃₃,K₁₂,K₂₃,K₁₃)^T

b＝(K_g(n₁),K_g(n₂),...,K_g(n_n))^T

式中：n_n1,n_n2,n_n3为渗透系数方向单位矢量在x,y,z轴上分量；

K₁₁,K₂₂,K₃₃,K₁₂,K₂₃,K₁₃为三维渗透张量的六个独立分量；

K_g(n_i)为第i个方向的渗透系数；

根据最小二乘法求解出渗透系数张量的六个分量x，获得三维渗透张量K：

$K=\left[\begin{array}{ccc}{K}_{11}& K_{12}& K_{13}\\ K_{12}& K_{22}& K_{23}\\ K_{13}& K_{23}& K_{33}\end{array}\right].$

2.根据权利要求1所述的一种三维渗透张量计算方法，其特征在于，所述的步骤二中，迹长和方向角所服从的密度分布参数包括迹长的均值和标准差，方向角的均值和标准差。

3.根据权利要求2所述的一种三维渗透张量计算方法，其特征在于，所述的步骤二中，二位裂隙生成域的尺寸为裂隙全局域尺寸的一半，裂隙全局域尺寸为裂隙迹长的均值的2-4倍。

4.根据权利要求1所述的一种三维渗透张量计算方法，其特征在于，所述的步骤二中，所述的视倾角γ转换通过tanγ＝tanαcosω计算，其中γ为视倾角，α为真倾角，ω为视倾向即剖面方向与倾向之夹角。