裂隙介质渗透性的升尺度转换研究.pdf

第4 2卷 第4期2 0 2 3年 7月 地 质 科 技 通 报B u l l e t i n o f G e o l o g i c a l S c i e n c e a n d T e c h n o l o g y V o l.4 2 N o.4J u l.2 0 2 3董晓飞,胡成,曹孟雄,等.裂隙介质渗透性的升尺度转换研究J.地质科技通报,2 0 2 3,4 2(4):2 5 9-2 6 7.D o n g X i a o f e i,H u C h e n g,C a o M e n g x i o n g,e t a l.S t u d y o n t h e u p s c a l i n g t r a n s f o r m a t i o n o f h y d r a u l i c c o n d u c t i v i t y i n f r a c t u r e d m e d i aJ.B u l l e t i n o f G e o l o g i c a l S c i e n c e a n d T e c h n o l o g y,2 0 2 3,4 2(4):2 5 9-2 6 7.基金项目:国家自然科学基金项目(4 2 0 2 2 0 1 8)作者简介:董晓飞(2 0 0 0),女,现正攻读水利工程专业硕士学位,主要从事裂隙渗流与数值模拟等方面的研究工作。E-m a i l:2 7 4 3 3 8 0 6 2 0c u g.e d u.c n通信作者:陈 刚(1 9 6 7),男,副教授,主要从事数值模拟及环境地质方面的研究工作。E-m a i l:c h e n g a n g c u g.e d u.c n E d i t o r i a l O f f i c e o f B u l l e t i n o f G e o l o g i c a l S c i e n c e a n d T e c h n o l o g y.T h i s i s a n o p e n a c c e s s a r t i c l e u n d e r t h e C C B Y-N C-N D l i c e n s e.裂隙介质渗透性的升尺度转换研究董晓飞,胡 成,曹孟雄,张 涛,陈 刚(中国地质大学(武汉)环境学院,武汉 4 3 0 0 7 8)摘 要:研究裂隙介质渗透性的升尺度转换,对准确刻画裂隙介质渗流场特征具有非常重要的意义。以某地下水封洞库的结晶岩裂隙介质统计数据为基础,应用M o n t e-C a r l o随机模拟技术生成二维离散裂隙网络(D F N)模型,计算得到变尺寸模拟域的渗透性参数及不同尺寸网格化划分后各网格单元的等效渗透系数,并对网格单元等效渗透系数进行升尺度运算。结果表明:模拟域尺寸达到渗透性典型单元体(R E V)尺寸2 2 m2 2 m后,模拟域可视为等效连续介质;网格化处理后,小于R E V尺寸的网格单元升尺度运算得到的等效渗透系数显著小于裂隙网络模型计算得到的对应复合网格单元的等效渗透系数。因此,渗流计算单元尺寸达到R E V尺寸后,其渗透性参数可以有效代表研究区内更大尺寸区域的渗透性特征;当渗流计算单元尺寸小于R E V尺寸时,其渗透性参数无法有效代表研究区内更大尺寸区域的渗透性特征,此时对渗透性参数进行参数升尺度运算往往具有低估的误差,不具有实际意义。关键词:离散裂隙网络模型;升尺度;R E V;渗透张量;裂隙介质2 0 2 3-0 1-1 4收稿;2 0 2 3-0 5-2 9修回;2 0 2 3-0 5-3 0接受中图分类号:X 1 4 1 文章编号:2 0 9 6-8 5 2 3(2 0 2 3)0 4-0 2 5 9-0 9d o i:1 0.1 9 5 0 9/j.c n k i.d z k q.t b 2 0 2 3 0 0 2 3 开放科学(资源服务)标识码(O S I D):S t u d y o n t h e u p s c a l i n g t r a n s f o r m a t i o n o f h y d r a u l i c c o n d u c t i v i t y i n f r a c t u r e d m e d i aD o n g X i a o f e i,H u C h e n g,C a o M e n g x i o n g,Z h a n g T a o,C h e n G a n g(S c h o o l o f E n v i r o n m e n t S t u d i e s,C h i n a U n i v e r s i t y o f G e o s c i e n c e s(Wu h a n),Wu h a n 4 3 0 0 7 8,C h i n a)A b s t r a c t:O b j e c t i v eI t i s v e r y i m p o r t a n t t o s t u d y t h e u p s c a l i n g t r a n s f o r m a t i o n o f h y d r a u l i c c o n d u c t i v i t y i n f r a c t u r e m e d i a f o r a c c u r a t e l y c h a r a c t e r i z i n g t h e s e e p a g e f i e l d c h a r a c t e r i s t i c s.M e t h o d sB a s e d o n t h e f r a c-t u r e m e d i u m s t a t i s t i c s o f c r y s t a l l i n e r o c k i n a n u n d e r g r o u n d w a t e r-s e a l e d c a v e r n,a 2-d i m e n s i o n a l d i s c r e t e f r a c t u r e n e t w o r k(D F N)m o d e l i s g e n e r a t e d u s i n g t h e M o n t e-C a r l o s t o c h a s t i c s i m u l a t i o n t e c h n i q u e.T h e h y d r a u l i c c o n d u c t i v i t y p a r a m e t e r s o f t h e v a r i a b l e s i z e s i m u l a t i o n d o m a i n a n d t h e e q u i v a l e n t h y d r a u l i c c o n-d u c t i v i t y o f e a c h g r i d e l e m e n t a f t e r m e s h i n g i n d i f f e r e n t s i z e s a r e c a l c u l a t e d.T h e h y d r a u l i c c o n d u c t i v i t y r e p r e s e n t a t i v e e l e m e n t a r y v o l u m e(R E V)o f t h e s t u d y a r e a w a s a n a l y s e d b y t h e v a r i a t i o n i n h y d r a u l i c c o n-d u c t i v i t y p a r a m e t e r s w i t h t h e s i z e o f t h e s i m u l a t i o n a r e a,a n d t h e e q u i v a l e n t h y d r a u l i c c o n d u c t i v i t y o f g r i d c e l l s s m a l l e r t h a n R E V w a s c a l c u l a t e d b y u p s c a l i n g.R e s u l t sT h e r e s u l t s s h o w t h a t t h e s i m u l a t i o n d o m a i n c a n b e r e g a r d e d a s a n e q u i v a l e n t c o n t i n u u m w h e n t h e s i z e o f t h e R E V r e a c h e s 2 2 m2 2 m.A f t e r m e s h i n g t r e a t m e n t,t h e e q u i v a l e n t h y d r a u l i c c o n d u c t i v i t y o f t h e m e s h c e l l s w h i c h i s s m a l l e r t h a n R E V c a l c u l a t e d b y h t t p s:/d z k j q b.c u g.e d u.c n 地质科技通报 2 0 2 3年 t h e u p s c a l i n g o p e r a t i o n i s s i g n i f i c a n t l y s m a l l e r t h a n t h a t o f t h e c o r r e s p o n d i n g c o m p o s i t e m e s h c e l l s c a l c u-l a t e d b y t h e c r a c k n e t w o r k m o d e l.C o n c l u s i o nT h e r e f o r e,w h e n t h e s i z e o f t h e s e e p a g e c a l c u l a t i o n u n i t r e a c h e s t h e R E V s i z e,i t s h y d r a u l i c c o n d u c t i v i t y p a r a m e t e r s c a n e f f e c t i v e l y r e p r e s e n t t h e h y d r a u l i c c o n d u c-t i v i t y c h a r a c t e r i s t i c s o f a l a r g e r a r e a i n t h e s t u d y a r e a.H o w e v e r,w h e n t h e s i z e o f t h e s e e p a g e c a l c u l a t i o n u n i t i s s m a l l e r t h a n t h e R E V s i z e,i t s h y d r a u l i c c o n d u c t i v i t y p a r a m e t e r s c a n n o t e f f e c t i v e l y r e p r e s e n t t h e h y d r a u l i c c o n d u c t i v i t y c h a r a c t e r i s t i c s o f t h e l a r g e r a r e a i n t h e s t u d y a r e a.A s a c o n s e q u e n c e,p a r a m e t e r u p-s c a l i n g c a l c u l a t i o n o n t h e h y d r a u l i c c o n d u c t i v i t y p a r a m e t e r s o f t e n h a s u n d e r e s t i m a t e d e r r o r s a n d i s n o t o f p r a c t i c a l s i g n i f i c a n c e.F o r f r a c t u r e d m e d i a s t u d y a r e a s w i t h i n s u f f i c i e n t d a t a,i t i s o f t e n d i f f i c u l t t o d e t e r-m i n e t h e h y d r a u l i c c o n d u c t i v i t y R E V o f t h e s t u d y a r e a.I n t h i s c a s e,i t c a n b e c o n s i d e r e d t h a t t h e r e i s u s u-a l l y a n u n d e r e s t i m a t i o n e r r o r w h e n t h e h y d r a u l i c c o n d u c t i v i t y o f t h e s m a l l-s c a l e a r e a o b t a i n e d f r o m t h e h y d r o g e o l o g i c a l t e s t i s u p s c a l e d.T h i s c o n c l u s i o n p r o v i d e s a t h e o r e t i c a l b a s i s f o r t h e n u m e r i c a l s i m u l a t i o n o f t h e s e e p a g e f i e l d i n f r a c t u r e m e d i a i n v a r i o u s r e l a t e d e n g i n e e r i n g p r o j e c t s.K e y w o r d s:d i s c r e t e f r a c t u r e n e t w o r k m o d e l;u p s c a l i n g;R E V;h y d r a u l i c c o n d u c t i v i t y t e n s o r;f r a c t u r e d m e d i aR e c e i v e d:2 0 2 3-0 1-1 4;R e v i s e d:2 0 2 3-0 5-2 9;A c c e p t e d:2 0 2 3-0 5-3 0 尺度效应最初是在研究渗透系数和溶质运移弥散度的过程中发现的一种现象:渗透系数和弥散度的实验室测量值往往比野外实测值小好几个数量级,这种渗透系数等参数随着研究尺度的变化而变化的现象称为尺度效应1-2。一般而言,空间上分布的参数在高度非均质的条件下都具有强烈的尺度效应。目前,国内外对渗透性尺度效应的研究手段主要分为物理试验法和数值试验法。孙蓉琳等3在溪洛渡坝区的玄武岩区域进行了多组不同尺度的水文地质试验,发现大尺度试验下得到的玄武岩渗透系数比小尺度试验更大且数值分布更为集中。刘日成等4基于UD E C平台进行二次开发,研究了裂隙网络模型的尺度效应并确定了典型单元体(R E V)尺度。高超等5基于C OM S O L数值模拟软件研究了单裂隙渗透性的尺度效应,结果表明单裂隙的渗透性随试样尺寸增大而增大,达到某一临界值后稳定。在裂隙介质渗流场模拟过程中,由于尺度效应的存在,当试验测量尺度显著小于模拟尺度时,试验测量尺度得到的渗透性参数值往往不能直接运用于模拟尺度。因此,升尺度方法被提出来用以解决渗透性参数尺度转换的问题,其目的在于将精细尺度(试验尺度)下的参数转换为粗尺度(模拟尺度)研究中的参数,并同时保持该参数场的有效性。升尺度方法通常包括统计平均法、体积平均法、均匀化理论和重整化技术,而根据对象的不同,又可分为参数升尺度和过程升尺度2种:参数升尺度指对于地下水中溶质的有效迁移和反应参数升尺度;过程升尺度指迁移方程升尺度6-9。覃荣高等1 0对升尺度转换问题进行系统综述,指出了升尺度转换在目前水文地质研究中主要解决的问题以及存在的问题。郭富欣等1 1对不同升尺度方法运算前后的模型分别进行数值模拟,通过数值模拟结果对比进行了升尺度方法优选以确定最适合研究区的升尺度方法。由于结晶岩裂隙介质往往具有高度的非均质性和空间随机性,且岩石基质自身渗透性较低,其渗透性的空间尺度效应表现得尤为显著1 2。因此对结晶岩裂隙介质渗透性升尺度转换的研究,具有十分重要的意义。然而,针对参数升尺度,大部分学者提出的升尺度算法是针对非均质孔隙介质的,裂隙介质的升尺度算法少,且大多是基于一定的研究场地,算法非常复杂,并不具有普适性,而升尺度的目的正是避免繁琐的计算,所以这些算法对于大型工程项目的实用性比较有限1 3-1 5。针对目前存在的一些问题和实际工程需求,本研究将应用普适性的升尺度方法对结晶岩裂隙介质渗透性的升尺度转换做更为详细的探讨。笔者拟以某地下水封洞库工程为例,在前期高精度勘察工作的基础上分析研究区裂隙发育规律,采用M o n t e-C a r l o方法建立研究区岩体二维离散裂隙网络模型。通过编译对应MAT L A B程序分别计算不同尺寸模拟域的渗透性参数,研究其随模拟域尺寸的变化规律并分析渗透性R E V尺寸,之后对裂隙网络模型进行不同尺寸网格化处理,计算各个网格单元的等效渗透系数,并采用3种普适性的升尺度方法对小于R E V尺寸的网格单元的等效渗透系数进行运算,并将计算结果与裂隙网络模型计算得到的对应复合网格单元的等效渗透系数进行对比,从而为实际模型应用中渗流场模拟提供借鉴。1 研究区二维离散裂隙网络模型1.1研究区概况研究区地下水封洞库工程建设场地位于山东省烟台市临港经济技术开发区内,在大地构造上处于062 第4期董晓飞等:裂隙介质渗透性的升尺度转换研究华北断块区东部的胶辽断块之中,地层岩性划分成3个种类:第四系冲洪积层(Q h l,Q h b,Qs);燕山早期中粗粒黑云二长花岗岩(2(1)5);煌斑岩脉、石英脉等各种岩脉。库区选址已避开周边主要断裂,但洞库建设区域内仍发育多条次级断裂,破碎带以及节理裂隙密集带。1.2二维离散裂隙网络模型生成基于M o n t e-C a r l o法对二维离散裂隙网络进行模拟的前提是单个裂隙能够用裂隙的产状、迹长、隙宽以及中心点坐标等几何参数完全表征,即在确认这些几何参数以后,可以在裂隙网络模型中完全刻画出这条裂隙的所有信息1 6-1 9。以研究区地下水封洞库的结晶岩裂隙介质统计数据为基础,主要是通过地表裂隙测量、水幕巷道地质素描图以及钻孔电视智能成像等数据资料对研究区内的裂隙几何特征参数进行统计和分析。研究区的地下裂隙岩体共发育3组较为明显的优势结构面,第1组裂隙结构面的倾向为N E 0 8 0,倾角为3 0 4 5;第2组裂隙结构面的倾向为S E 1 0 0 1 4 0,倾角为3 5 5 0;第3组裂隙结构面的倾向为NW 2 9 0 3 5 0,倾角为7 0 8 5。3个优势裂隙组的裂隙方向角都近似服从正态分布(图1),迹长都近似服从对数正态分布(图2)。由于岩体裂隙宽度大部分都在0.0 0 1 m以内,现有的技术手段都很难对这些裂隙宽度进行精确的测量,因此本研究参考前人给出的经验值,假定隙宽服从均值为0.0 0 0 4 m、标准差为0.0 0 0 2 m的对数正态分布。3组裂隙的方向角、迹长、隙宽和密度的统计规律如表1所示。此外,由于裂隙的分布具有高度的空间随机性,结合地质素描图中所观察到的裂隙发育情况,假设在二维裂隙网络模拟中,裂隙中心点的位置服从均匀分布。图1 优势结构面1,2,3方向角分布直方图F i g.1 H i s t o g r a m o f t h e a n g u l a r d i s t r i b u t i o n i n d i r e c t i o n s 1,2 a n d 3 o f t h e d o m i n a n t s t r u c t u r a l p l a n e图2 优势结构面1,2,3迹长分布直方图F i g.2 H i s t o g r a m o f t r a c e l e n g t h d i s t r i b u t i o n o f d o m i n a n t s t r u c t u r a l p l a n e s 1,2 a n d 3表1 优势裂隙组各个几何参数分布类型及统计规律T a b l e 1 D i s t r i b u t i o n t y p e s a n d s t a t i s t i c a l r u l e s o f g e o m e t r i c p a r a m e t e r s i n t h e d o m i n a n t f r a c t u r e g r o u p优势裂隙组数方向角(正态分布)迹长(对数正态分布)隙宽(对数正态分布)均值/()标准差/()均值/m标准差/m均值/m标准差/m密度/(条m-2)11 4 2.2 93.7 85.2 32.9 00.0 0 0 40.0 0 0 20.2 924 5.9 15.5 15.6 22.6 70.0 0 0 40.0 0 0 20.2 338 1.2 93.6 36.4 23.3 40.0 0 0 40.0 0 0 20.2 4162h t t p s:/d z k j q b.c u g.e d u.c n 地质科技通报 2 0 2 3年 根据表1中的迹长数据,迹长均值为5.2 36.4 2 m,为确保后续裂隙渗流计算数据稳定,确定裂隙生成域为6 0 m6 0 m,计算域取3 0 m3 0 m。在MAT L A B环境中编写相应的程序,并将表1中裂隙各个几何参数的概率分布模型参数代入,生成符合实际场地裂隙发育规律的二维裂隙网络模型,如图3。其中第1组优势裂隙为1 0 4 4条,第2组优势裂隙为8 6 4条,第3组优势裂隙为8 2 8条,生成域内共2 7 3 6条裂隙,与表1优势裂隙组密度分布情况一致。图3 基于烟台洞库数据的二维裂隙网络模型效果图F i g.3 E f f e c t o f t h e t w o-d i m e n s i o n a l f r a c t u r e n e t w o r k m o d e l b a s e d o n Y a n t a i c a v e d a t a2 裂隙岩体渗透张量确定离散裂隙网络渗流模型忽略了低渗透性的岩块孔隙系统,假定流体只在裂隙网络中运动,从而将裂隙岩体介质系统转化为裂隙网络介质系统。在左、右边界为存在水头差的定水头边界,水头值分别为1 5,1 0 m,上下边界为隔水边界的条件下,本研究基于离散裂隙网络稳定渗流数学模型计算裂隙岩体渗透系数,并通过求解不同方向上的渗透系数来拟合渗透椭圆的方法,得到裂隙岩体的渗透张量。左右定水头边界上的渗流量公式为:Q=-A3T AT1H1-A3T AT2H2-A3T AT3H3(1)式中:Q为渗流量;A1,A2,A3为衔接矩阵,分别描述裂隙网络内节点、裂隙与上下边界交点和裂隙与左右边界交点与线单元的衔接关系;T为对角矩阵,该矩阵中对角线上的元素为Tj=g b3jhjLj,其中,为裂隙密度,g为重力加速度,是与裂隙面粗糙程度相关的系数,本研究假设裂隙光直平滑,取1/1 2,bj为第j条裂隙的隙宽,为地下水的动力黏滞系数;hj为水头,Lj为渗透途径。H1,H2,H3为水头矩阵,令渗流域内的裂隙网络内节点数为N1,第二类边界点数量为N2,第一类边界点数量为N3,则H1,H2,H3分别为N11,N21,N31阶水头向量。根据达西定律,沿着水头梯度方向上的渗透系数K可由下式求得:K=vJ=Q/(L1)(H1-H2)/L=QH1-H2(2)式中:K为渗透系数;v为渗透流速;J为水力梯度;H1、H2为上下游过水断面的水头。对于二维裂隙网络模型渗流分析问题,渗透张量Ki j是一个二阶对称张量。将结构面的渗透性平均到裂隙岩体中,得到均质各向异性的裂隙介质,Ki j为渗透张量,记作K。用矩阵可表示为:K=Kx xKx yKy xKy y(3)图4 裂隙网络模型绕中心点旋转示意图F i g.4 D i a g r a m o f t h e f r a c t u r e n e t w o r k m o d e l r o t a t i n g a-r o u n d t h e c e n t e r p o i n t对于本研究的二维裂隙网络稳定渗流问题,求解裂隙介质渗透张量Ki j的具体过程如下:根据上文所述的离散裂隙网络稳定渗流数学模型方法,在给定的边界条件下计算出计算域沿水力梯度方向上的等效渗透系数。绕模型中心点旋转裂隙网络,设置为每隔3 0 旋转一次,共取得1 2个不同方向的计算域网络模型,并分别计算其沿水力梯度方向上的渗透系数,见图4。以1 2个方向的度数及对应角度所计算出的渗透系数的平方根倒数为坐标点,在极坐标系下采用最小二乘法将这1 2个点拟合成渗透椭圆,见图5。求出渗透椭圆的长半轴和短半轴长度,从而得到裂隙介质的渗透主值K1、K2,并根据渗透主值与极坐标的0 方向夹角确定渗透主方向,然后由下式计算出渗透张量的分量值:Kx x=12(K1+K2)+12(K1-K2)c o s 2(4)262 第4期董晓飞等:裂隙介质渗透性的升尺度转换研究 Ky y=12(K1+K2)-12(K1-K2)c o s 2(5)Kx y=Ky x=12(K1-K2)s i n 2(6)图5 裂隙网络模型渗透椭圆示意图F i g.5 H y d r a u l i c c o n d u c t i v i t y e l l i p s e d i a g r a m o f t h e f r a c-t u r e n e t w o r k m o d e l为了判断计算出的渗透张量能否较准确地表征裂隙网络模型的渗流特征,需要进行均方根误差RMS的计算。RMS是评价拟合过程中渗透系数计算值和渗透椭圆拟合值的拟合程度的一个重要指标,RMS越小代表拟合程度越好,其计算方法如下:RMS=2K1+K2 1NN1Kg()-Kc()2(7)式中:N为不同方向计算出的渗透系数个数,取1 2;Kg()为方向角为时的渗透系数计算值;Kc()为方向角为时应用最小二乘法拟合出的渗透系数。3 渗透性R E V分析R E V是能够从宏观角度反映出研究区某种属性的最小区域。对于本次的裂隙网络稳定渗流模拟研究来说,R E V指的是渗透性典型单元体。本研究将根据变尺寸模拟域的渗透性相关参数:渗透张量分量Kx x,Ky y、均方根误差RMS、渗透张量分量的变异系数C Vx x,C Vy y随模拟域尺寸变化情况(如图68所示)分析研究区渗透性R E V。其中,为避免分析结果受到M o n t e-C a r l o随机模拟技术生成的裂隙网络模型具有一定的随机性的影响,基于裂隙几何参数统计数据,重复生成了1 0组6 0 m6 0 m的裂隙网络模型(D F N 1D F N 1 0),计算域取为3 0 m3 0 m,分别计算了1 0组D F N模型在变尺寸模拟域(4 m4 m,6 m6 m,3 0 m3 0 m)中的渗透性参数并进行分析。图6 渗透张量分量随模拟域尺寸变化图F i g.6 H y d r a u l i c c o n d u c t i v i t y t e n s o r c o m p o n e n t c h a n g e s w i t h t h e s i z e o f t h e s i m u l a t i o n d o m a i n图7 RMS随模拟域尺寸变化图F i g.7 RMS v a r i a t i o n w i t h t h e s i z e o f t h e s i m u l a t i o n d o m a i n 当Kx x、Ky y的数值基本稳定时,该裂隙网络模型的渗透性不再发生明显变化,认为此时研究尺寸达到了R E V尺寸;根据 h m a n等2 0的研究,当RMS0.2时,可认为此时的研究尺寸达到了R E V尺寸;对于C Vx x,C Vy y,设定一个容许值,当变异系数小于容许值时,可认为此时的研究尺寸达到了R E V尺寸,根据M i n等2 1的研究,变异系数的容许值可取为1 0%。结合图68对渗透性参数Kx x,Ky y,RMS,C Vx x,C Vy y随着模拟域尺寸增加而变化的规律进行分析,可以发现1 0组离散裂隙网络模型的渗透性362h t t p s:/d z k j q b.c u g.e d u.c n 地质科技通报 2 0 2 3年 图8 变异系数随模拟域尺寸变化图F i g.8 C o e f f i c i e n t o f v a r i a t i o n w i t h t h e s i z e o f t h e s i m u l a-t i o n d o m a i n参数计算结果均显示出了相似的规律。在模拟域尺寸较小时,不同随机裂隙网络计算出的渗透性参数差异较大,各向异性显著;随着模拟域尺寸的增大,渗透性参数差异变小,且都在模拟域边长达到2 0 m或2 2 m时达到R E V判定临界值,因此考虑到裂隙网络模型生成的随机性,可将该研究区的R E V尺寸取为2 2 m2 2 m。在达到R E V尺寸后,模拟域可等效为连续介质处理,模拟域计算得到的渗透性参数可以代表研究区的渗透性特征。因此,渗流计算单元尺寸大于R E V尺寸时,渗透性参数趋于一致,可以直接应用R E V渗透性参数作为研究区内更大尺寸区域的等效渗透性参数进行裂隙介质渗流场模拟。4 等效渗透系数升尺度转换本研究进行升尺度运算时没有考虑渗透主值与渗透主方向的影响,所以升尺度运算仅对裂隙网络旋转角度为0 时的等效渗透系数进行。由于升尺度运算采用4个网格单元的等效渗透系数进行最为合适,所以需改变裂隙网络模型生成域尺寸便于进行升尺度运算。基于裂隙几何参数分布规律及统计数据,设定裂隙网络模型生成域为8 0 m8 0 m,计算域为4 0 m4 0 m。网格单元尺寸极小时,裂隙网络连通性差,等效渗透系数计算结果没有规律性,故取最小网格单元尺寸为5 m5 m。对裂隙网络模型进行不同尺寸网格化处理(5 m5 m、1 0 m1 0 m、2 0 m2 0 m),并对5 m5 m、1 0 m1 0 m网格单元的等效渗透系数进行升尺度运算。为了便于区分,将升尺度运算得到的等效渗透系数记为K,5 m5 m网格单元升尺度运算得到的等效渗透系数记为K5,1 0 m1 0 m网格单元升尺度运算得到的等效渗透系数记为K1 0。裂隙网络模型计算得到的1 0 m1 0 m复合网格单元的等效渗透系数记为K1 0,2 0 m2 0 m复合网格单元的等效渗透系数记为K2 0。4 0 m4 0 m的计算域可划分为6 4个5 m5 m的网格单元,1 6个1 0 m1 0 m的网格单元,4个2 0 m2 0 m的网格单元,分别计算出1 0组裂隙网络模型不同尺寸网格单元的等效渗透系数,采用4个5 m5 m网格单元的等效渗透系数进行升尺度运算以得到K5与K1 0对应,4个1 0 m1 0 m网格单元的等效渗透系数进行升尺度运算以得到K1 0与K2 0对应。本研究采用参数升尺度方法中的3种普适性的升尺度方法对5 m5 m网格单元与1 0 m1 0 m网格单元的等效渗透系数进行升尺度运算,这3种升尺度方法无需确定流动状态,可直观展现出网格单元的等效渗透系数是否具有代表性。4.1统计平均法一个已知渗透系数分布的区块的等效渗透系数是有上限和下限的,其区间的范围由 W i e n e r 界限确定(调和平均数=下限,算术平均数=上限)2 2。通常使用加权平均来描述 W i e n e r界限之间的等效渗透系数范围,见式(8)。Ke f f=1/p=1VVK(x)pdV 1/p(8)式中:Ke f f为等效渗透系数;V为区块体积;p为幂指数,指数范围为-1到1。当p=-1时,上式表示调和平均数Kh(下限),当p=1时,上式表示算数平均数Ka(上限),当p无限趋近于0时,上式表示几何平均数Kg。4.2标准重整化我们称K i n g2 3提出的升尺度方法为标准重整化,通过类比介质中流体流动与通过电阻的电流流动而提出的一种二维升尺度方法。二维标准重整化升尺度的等效渗透系数Ks是通过一系列连续的4个单元聚合来计算的,如图9。通过下式来对二维网格单元的等效渗透系数进行升尺度运算:Ks=4 K1+K2+K3+K4(K1+K3)(K2+K4)+3(K1+K2)(K3+K4)K2K4(K1+K3)+K1K3(K2+K4)(9)4.3简化重整化简化重整化方法是在综合了C a r d w e l l&P a r-s o n s边界及重整化方法原理的基础上通过迭代算法对相邻单元的等效渗透系数进行组合,算法具体步骤为1 1,2 4:在第一步中,根据计算方向交替地以并联和串联的方式对单元进行分组。如果2个单元462 第4期董晓飞等:裂隙介质渗透性的升尺度转换研究图9 二维标准重整化过程示意图F i g.9 S c h e m a t i c d i a g r a m o f a t w o-d i m e n s i o n a l s t a n d-a r d r e n o r m a l i z a t i o n p r o c e s s是串联的,它们被1个唯一的单元取代,等效渗透系数是单元渗透系数的调和平均值Kxh。如果2个单元是并联的,等效渗透系数是单元渗透系数的算术平均值Kya。这个基本过程系统地重复,直到得到1个唯一的值。分组的顺序影响最终结果。在二维中,可以从沿着x方向的串联分组开始,然后沿着y方向并联分组,重复这个基本算法,直到得到一个表示Cm i n的值,见式(1 0);或者可以从沿着y方向并联分组开始,然后沿着x方向串联分组,重复这个算法以得到一个表示Cm a x的值,见式(1 1)。具体二维简化重整化过程如图1 0所示。Cm i n=Kya(Kya(Kxh)(1 0)Cm a x=Kxh(Kxh(Kya)(1 1)图1 0 二维简化重整化过程示意图F i g.1 0 S c h e m a t i c d i a g r a m o f a t w o-d i m e n s i o n a l s i m p l i f i e d r e n o r m a l i z a t i o n p r o c e s s5 结果与讨论分别计算出1 0组裂隙网络模型不同尺寸网格单元(5 m5 m,1 0 m1 0 m,2 0 m2 0 m)的等效渗透系数,通过以上3种升尺度方法对1 0组裂隙网络模型5 m5 m网格单元与1 0 m1 0 m网格单元的等效渗透系数