热流体数值模拟模型

Question

如题所述

Answer 1

(1)数学模型

根据上述非均质、水平方向各向同性、垂直方向存在变异的空间三维结构、非稳定热流体系统,可用如下微分方程的定解问题来描述:

图5-16 Ng组三维实体模型建模演示图

沉积盆地型地热田勘查开发与利用

式中:Ω为渗流区域;h为热流体的水位标高(m);K_x,K_y,K_z分别为热流体储积层中x,y,z方向的渗透系数(m/d);K_n为边界面法向方向热流体储积层渗透系数(m/d);K_z→为边界内侧热流体储积层垂向渗透系数(m/d);S为储水率(1/m);σ为通用水头边界的阻力系数,σ=L/K₁,L为模型边界到通用水头边界的水平距离(m),K₁为模型边界到通用水头边界之间热流体储积层平均渗透系数(m/d);ε为源汇项(1/d);h₀为初始水位分布(m),h₀=h₀(x,y,z);h_b为通用水头边界上的水位分布(m);

为边界面的法线方向;

为热储层顶、底边界面的法线方向;Γ₁为通用水头边界;Γ₂为渗流区域的侧向隔水边界;Γ₃为渗流区域顶、底边界;

(2)数值模型

绝大部分数学模型是无法用解析法求解的。数值化就是将数学模型转化为可解的数值模型。常用数值化有限单元法和有限差分法。

本次工作根据热储层地热地质概念模型特征、实际资料丰富程度以及精度要求,应用有限差分法对馆陶组热储进行评价。

模型软件选用模拟地下水流动和污染物运移等特性的计算机程序——Modflow (The modular finite-difference groundwater flow model)。该系统是三维地下水流和溶质运移模拟评价软件。

(3)计算区单元剖分

在进行热流体动态模拟前,首先必须对模拟区进行单元剖分,剖分时除了遵循一般的剖分原则外,还应充分考虑如下实际情况:①计算区的边界、岩性分区边界、行政分区边界等;②观测孔尽量放在剖分单元的结点上;③集中开采区域,由于水力坡度较大,流场变化趋势较大,所以剖分时适当加密。因此在平面上采用不等距矩形剖分法,在开采井集中区域适当加密剖分间距。沿边界将外部和缺失区域定义为非活动单元格。剖分为60行、55列,其中有1580个非活动单元。最大剖分单元约1.8×1.8km²正方形单元,最小剖分单元为0.28×0.38km²长方形单元。而在垂向上热储层按NgⅠ、中间弱透水层NgⅡ和NgⅢ进行剖分。则Ng有(60×55-1580)×3三层一共是5160个剖分单元,具体情况见图5-17和图5-18。

(4)边界条件

Ng热储层可以概化出3种类型的边界条件:即①自然流量边界;②零流量隔水边界;③沿沧东断裂形成的内部补给边界(图5-19)。现对不同类型边界条件的处理方法分述如下:

A.自然流量水头边界

自然流量边界是利用MODFLOW软件中的变水头边界程序包(GHB)中的通用水头边界来处理。使用这个边界条件的目的是避免不必要的向外扩展模型区域,以满足要素影响模型中的水头。因此,变水头边界条件通常被指定在模型域的外面。从外部源进入或者流出计算单元的热流体量与该计算单元水头和外部源水头之差成正比。

在通用水头边界单元上,从单元外部的热流体进入或流出计算单元的热流体量与该计算单元水头h和通用水头边界热流体水头h_b之差成正比。据此,可确立计算单元热流体量与计算单元水头间的线性关系,即:

图5-17 Ng平面单元剖分图

图5-18 Ng3D单元剖分图

图5-19 Ng计算边界条件图

沉积盆地型地热田勘查开发与利用

式中:Q_b为从通用水头边界进入计算单元的流量(m³/d);C_b为通用水头边界与计算单元间的水力传导系数(m²/d);(L×W)为边界与网格单元的截面积(m²);K为从模型网格中分离源/汇的含水层物质的平均渗透系数(m/d);D为边界源/汇到模型网格间的距离(m);h_b为通用水头边界上的热流体水头(m);h为计算单元热流体水头(m)。

当C_b很大时,GHB子程序包的作用与定水头边界相同。因此通用水头边界上的水位和水力传导系数C_b决定了边界的流入流出量Q_b的大小。

B.零流量隔水边界

除了新近系馆陶组西侧边界和部分缺失区域定为零流量隔水边界,各热储层顶底板也定为隔水边界。毫无疑问,对于隔水边界,其计算单元与外部域的流入流出量Q_f值为零。

C.沿沧东断裂形成的内部补给边界

由沧东断裂特性,沿沧东断裂形成了连通上下各热储层地热流体的一个通道。滨海塘沽区、大港区对Ng热流体的过量开采可导致该层热流体压力水头的过量下降,从而激发其下部各热储层对该层的补给,因此,沿沧东断裂一带应考虑其他含水层的垂向补给。本次工作通过数值模型软件 MODFLOW提供的河流子程序包来模拟沧东断裂对计算区热流体的补给量。其计算公式如下:

沉积盆地型地热田勘查开发与利用

式中:Q_fau为从断裂进入计算单元的流量(m³/d),进入取正值,流出取负值;C_fau为断裂带外部源与断裂带计算单元间的水力传导系数(m²/d);L、W为通过计算单元网格的断层长度和宽度(m);K为断裂破碎带间计算单元与下层热储层间平均渗透系数(m/d);M为断裂破碎带间计算单元与下层热储层间厚度(m);h_fau为断裂破碎带间计算单元下层热储层热流体水头(m);h为断裂破碎带间计算单元水头(m)。

沧东断裂带下层基岩热流体储积层水头一般都高于Ng热流体水头,断裂带的渗流补给量一般是正值。即下层基岩热流体通过断裂带顶托补给上层Ng热储层。但由于受原始资料限制,所求的流入流出量Q_fau精度较粗,不能准确给出。只能在流场模拟中,根据流场模拟结果来进行校准。

(5)模拟期的选择和时间离散

数值模拟自2002年1月始,至2006年10月结束,共5年10个月时间作为该层模型识别整个模拟时间长度。根据地热流体的实际开采利用特点,将一个开采周期年分为两个应力期(时段):即取暖期(11月10日至次年3月25日,共135天)和非取暖期(3月26日至11月9日,共230天),每个应力期又分成基本以15天为单位的时间步长。

(6)补给、径流、排泄条件

对于地热流体,由于埋藏较深,相应补排项比较单一:补给来源主要包括侧向径流补给、随着热储层热流体压力变化引起的热储层岩石骨架释水、部分同组层间越流补给量以及人工回灌量。新近系馆陶组过量开采可能通过沧东断裂激发部分基岩补给量;排泄方式主要为下游单元的侧向径流排泄、人工开采以及越流补给同组其他超采层。

对于侧向径流以及内部补给边界的处理前文已有介绍,在此不再赘述。以下主要分析其他补排项的处理。

A.弹性释水

引起岩石骨架释水的主要原因是热储层热流体压力变化,而岩石骨架释水量的大小,不但受压力变化的影响,热储层弹性释水系数的大小也是一个至关重要的参数。其值由以下公式计算而得

Q_i=L×W×S^*×Δh 5-25

式中:Q_i为单元弹性释水量(m³);L为单元长度(m);W为单元宽度(m);S^*为单元所在地热流体富集段弹性释水系数(无量纲);Δh为单元热流体压力水头变化量(m)

B.同组层间越流补给、排泄量

由于采用三维水流模型,层间越流量由MODFLOW模型软件根据所赋相关参数自动进行处理。

C.泥岩释水

泥岩释水是指土体中无所不在的结合水和毛细水由于土体压缩而释放的水量,泥岩释水可引起地面沉降。很明显,由于地热开采引起的地面沉降有限,泥岩释水量也不会太大,但也不应对其忽视。

据1998年12月《天津市塘沽区地热流体回灌研究》(李明朗,1998)指出,热流体开采量中黏性土释水量约占总释水量(弹性释水量+粘土释水)的1/5;而总释水量高时约占到热流体开采量40%左右。但随着资源开采的延续,黏性土层释水后压密度增加,其释水量应呈逐步衰减之势。而随着热储层压力下降,总释水量也会逐步减少。

D.人工开采、回灌量的处理

对于人工开采、回灌量,在模型中都是以开采井的形式给入,开采量为正值而回灌量为负值。

(7)水文地质参数

水文地质参数主要是指渗透系数K和弹性释水系数S^*,前者K体现了热储层流体的渗透能力;而后者S^*则体现了热储层储水(释水)能力。水文地质参数选用是依据三维模型的计算需要,主要采取分区赋值的方法。根据Ng热储层的埋藏分布和控热构造分布特征,利用所收集的降压试验资料进行综合统计分析,按所求水文地质参数进行初始分区(见图5-20)。有条件的区域项目前期补充地热井降压试验工作,空白区域根据邻区资料结合地层构造变化特征给出估计值。