一种三维渗透张量计算方法
【技术领域】
[0001] 本发明涉及一种三维渗透张量计算方法。
【背景技术】
[0002] 经典渗流理论是以连续介质假定为基础的，而实际情况中，大多数介质都是离散 的。为了能将离散介质适用于经典渗流理论，引入了表征单元体（REV)的概念，从而把实际 上离散的多孔介质等效成由表征单元体组成的连续介质。当表征单元体尺寸与研究区域尺 寸相比很小时，就可以用连续介质渗流理论进行分析。对于土体等孔隙介质，其表征单元体 非常小，可以直接运用连续介质渗流理论。对于岩体，由于岩体中通常会有各种节理裂隙， 在流体流经岩体时，主要是由裂隙通过的。裂隙和孔隙主要有两种差别：1.孔隙在三维方 向上延伸，且各个方向上的差别不大，而裂隙主要在二维方向上延伸，其两个方向上的尺寸 比第三个方向上的尺寸大很多。2.孔隙的延伸尺寸比研究域小得多，而裂隙的延伸尺度可 以达到整个研究域。
[0003] -直以来，工程师在计算中把岩体看成和土体一样的多孔介质，应用经典渗流理 论进行分析。而1959年法国Malpasset拱坝和1967年Vajiont拱坝事故的发生，使人们 认识到岩体渗流不能和土体一样直接采用连续介质渗流理论分析。
[0004] 如今，随着裂隙性油气田的开发、高放射性核废料的深埋处理和对拱坝岩体的开 发，使得裂隙岩体渗流的研究逐渐成为了热点问题。目前对裂隙岩体进行渗流计算时，主要 有：等效连续介质模型、离散裂隙网络模型和双重介质模型。但是由于后两者都有难以解 决的问题而影响了模型的应用：在双重介质模型中，裂隙网络和连续介质岩块耦合模型中 水交换量公式的准确性直接影响着模型的计算准确性，但是该模型中水交换量难以准确确 定；离散裂隙网络模型认为岩块不能导水而忽略岩块的导水性，但实际情况中，岩块的非饱 和渗透系数和裂隙渗透系数相比，一般不能忽略，故当考虑裂隙岩体非饱和渗流时，选用离 散裂隙模型并不合适；相比之下，等效连续介质模型分析裂隙岩体渗流，在理论和实际中均 有非常多的经验和成熟和计算程序可供借鉴。同时，该模型只需知道裂隙的几何信息的统 计值，而对其准确位置并不敏感，因此，在实际工程问题中多使用等效连续介质模型。渗透 张量是等效连续介质渗流计算模型中的一个关键参数，它直接反应了裂隙岩体非均质各向 异性的特点。近年来，许多学者和专家在等效连续介质渗流计算中，就如何准确确定渗透张 量这一问题进行了大量研究分析。
[0005] 在三维情况下，用等效连续介质模型分析裂隙岩体的渗流时，需要使用一个3阶 的渗透张量来表达裂隙渗流的各向异性特性。许多学者采用裂隙网络渗流的方法来计算三 维渗透系数张量。Long J C S(1985)提出圆盘裂隙网络模型。Dershowitz (1987)提出多 边形裂隙网络模型。Cacas M C(1990)提出圆形管道模型。这些方法的实质是将裂隙及其 交线看成一个点，将裂隙中的渗流简化成水在这些点之间的直线运动。但是，采用这些方法 进行渗流分析时，由于三维情况下裂隙网络的复杂性，需要花费大量精力在裂隙生成、筛选 以及裂隙与裂隙或边界的联系这些方面上。当岩体内部裂隙密度高时，可能在20米的尺度 范围内就能包含上万条对渗流有作用的裂隙，这时，对裂隙的生成、筛选以及相关计算工作 是非常巨大的。
【发明内容】
[0006] 为了解决目前在三维情况下分析裂隙岩体的渗流计算工作量过大的技术问题，本 发明提供一种可有效提高有效裂隙的计算效率的三维渗透张量计算方法。
[0007] 为了实现上述技术目的，本发明的技术方案是：
[0008] -种三维渗透张量计算方法，包括以下步骤：
[0009] 步骤一：输入由统计资料确定的各组裂隙面服从密度分布参数；
[0010] 步骤二：在三维坐标中假定一个二维裂隙生成域，将优势裂隙面产状进行视倾角 转换，获得裂隙在该二维生成域中迹长和方向角所服从的密度分布参数；
[0011] 步骤三：在二维生成域中根据步骤二得到的迹长和方向角所服从的密度分布参 数，采用ment-carlo法生成裂隙，并在生成域中心位置假定分析域尺寸，分析域的尺寸不 大于生成域的尺寸；
[0012] 步骤四：在分析域的范围内，对裂隙进行筛选，通过裂隙网络渗流模型计算分析域 二维面上与渗透水流方向垂直的边界节点的流量Q和水头H，计算分析域二维平面上水流 方向上的渗透系数K :
[0013]
⑴
[0014] 其中Q和Λ H为与渗透水流方向垂直的边界节点的流量和水头差；V为流速，J为 水力梯度，L为分析域尺寸，H1, H2为分析域区域左右两边的水头；
[0015] 步骤五：在步骤三的二维生成域中，将分析域旋转角度β后重复步骤四，直到旋 转一周为止，获得各次旋转后的渗透系数匕（0);
[0016] 步骤六：将各个方向上的渗透系数经过处理后拟合成椭圆，计算椭圆拟合质量 RMSnorn:
[0017]
(2)
[0018] 式中：N代表步骤五中计算得到的渗透系数的个数，kg(i3)是由旋转角度β后通 过裂隙网络模型计算得到的渗透系数；k n(i3)是拟合椭圆上相应方向角β上的渗透系数；
[0019] 步骤七：如果RMSnmiS 0. 2则认为椭圆质量好，其表征单元体REV为分析域尺寸； 反之，增大分析域尺寸，直到满足条件，但分析域尺寸大于生成域尺寸时，认为REV不存在；
[0020] 步骤八：以任一穿过生成域的直线为定轴旋转二维裂隙生成域，直至旋转360°， 重复步骤二~七，新旋转角度后的二维生成域的初始分析域尺寸取各REV中的最大值；
[0021] 步骤九：根据前面计算出的各个方向的渗透系数及其方向，组成下式：
[0022] Ax = b (3)
[0023] 其中：
[0024] CN 105138738 A ~P 3/6 页
[0025]
[0026]
[0027] 式中：nnl，nn2, nn3为的渗透系数方向的单位矢量的在x，y，z轴上分量；
[0028] Kn，K22, K33, K12, K23, K13为三维渗透张量的六个独立分量；
[0029] Kg(Ii1)为第i个方向的渗透系数；
[0030] 根据最小二乘法求解出渗透系数张量的六个分量X，获得三维渗透张量K :
[0031]
[0032] 所述的一种三维渗透张量计算方法，所述的步骤二中，迹长和方向角所服从的密 度分布参数包括迹长的均值和标准差，方向角的均值和标准差。
[0033] 所述的一种三维渗透张量计算方法，所述的步骤二中，二位裂隙生成域的尺寸为 裂隙全局域尺寸的一半，裂隙全局域尺寸为裂隙迹长的均值的2-4倍。
[0034] 所述的一种三维渗透张量计算方法，所述的步骤二中，所述的视倾角γ转换通过 tan γ = tana cosco计算，其中γ为视倾角，α为真倾角，ω为视倾向即剖面方向与倾向 之夹角。
[0035] 本发明的技术效果在于，在二维裂隙网络渗流分析的基础上，沿一定轴旋转二维 剖面，计算不同旋转角度下，二维剖面各个方向的渗透系数，并绘制椭圆球体，最后用回归 分析近似求得三维的渗透系数张量裂隙。由于在二维裂隙网络渗流计算中，采用一维的线 单元模拟裂隙，有效地避免了三维裂隙网络模型中复杂的裂隙处理工作。在裂隙计算时， 因为求解方程的维度由三维降低到二维，所以可以大幅度地提高有效裂隙的计算数量和效 率。
【具体实施方式】
[0036] 本方法包括以下步骤：
[0037] 步骤一：输入由统计资料确定的各组裂隙面服从密度分布参数；
[0038] 步骤二：在三维坐标中假定一个二维裂隙生成域，该二位裂隙生成域的尺寸为裂 隙全局域尺寸的一半，裂隙全局域尺寸为裂隙迹长的均值的2-4倍，将优势裂隙面产状进 行视倾角转换，其中，真倾向与视倾向之间夹角在倾斜面上斜交走向线所引的任一直线均 为视倾斜线（如在任一斜交岩层走向的露头断面上所见），视倾斜线与其在水平参考面上 的投影线的夹角，叫视倾角。这样可获得裂隙在该二维生成域中迹长和方向角所服从的密 度分布参数，当密度函数服从正态分布时，相关参数包括迹长的均值和标准差，方向角的均 值和标准差；标准差和均值为统计分析学中的公知含义；
[0039] 步骤三：在二维生成域中根据步骤二得到的迹长和方向角所服从的密度分布参 数，采用ment-carlo法生成裂隙，并在生成域中心位置假定分析域尺寸，分析域的尺寸不 大于生成域的尺寸；
[0040] 步骤四：在分析域的范围内，对裂隙进行筛选，由于分析域是一个矩形，有两边与 水流方向相同，两边与水流方向垂直，这里计算与水流方向垂直的边界上的节点，通过裂隙 网络渗流模