饱和多孔介质流固耦合渗流的数学模型

A辑第18卷第4期 水动力学研究与进展 Ser.A,Vol.18,No.4

2003年7月

JOURNALOFHYDRODYNAMICS July,2003

文章编号:1000-4874(2003)04-0419-08

饱和多孔介质流固耦合渗流的

数学模型

*

李培超, 孔祥言, 卢德唐

(中国科技大学力学和机械工程系,安徽合肥230027)

摘 要: 将基于多孔介质的有效应力原理引入流固耦合渗流中,并根据平衡条件得出了应力场方程;充分分析流固耦合渗流的物理特性,建立起孔隙度和渗透率动态模型;依据流体力学连续性方程,考虑流固耦合情形下多孔介质骨架变形特性和流体的可压缩性,得到了孔隙流体的连续性方程,即渗流场方程。在以上诸方程的基础上辅助以定解条件,建立起了完备的饱和多孔介质流固耦合渗流的数学模型。本文还对模型进行了简化分析,并与经典一维固结理论作定性对比,分析表明,经典的固结理论可认为是本文模型的一种简化。进一步的算例求解和验证工作正在进行中。

关 键 词: 饱和多孔介质;流固耦合;有效应力原理;孔隙度和渗透率动态模型;固结方程中图分类号: O357.3 文献标识码:A

Mathematicalmodelingofflowinsaturatedporousmedia

onaccountoffluid-structurecouplingeffect

LIPe-ichao, KONGXiang-yan, LUDe-tang

(DepartmentofMechanicsandMechanicalEngineering,UniversityofScienceand

TechnologyofChina,Hefei230027,China)

Abstract: Principleofeffectivestressinporousmediaisintroducedintofluid-structurecouplingtheory,thenaccordingtotheequilibriumcondition,thestressequationisestablished.Onaccountofphysicalpropertiesofcoupling,thedynamicmodelsofporosityandpermeabilityareacquired.Meanwhile,inconsiderationofthedeformationoftheskeletonandthecompressibilityofthefluid,theflowequationisobtained.Thus,sel-fcontainedgoverningequationsforfluidstructurecouplinginporousmediaareconstructedbycombiningtheaforementionedequationsandtheconstraintconditions.Onedimensionalconsolidationequationisobtainedbysimplify-ingthemodel,anditiscomparedwiththetraditionalconsolidationequationbyqualitativeanalysis.Theresultsindicatethatthetra-ditionalconsolidationtheorycanberegardedasoneparticularcaseofourmodel.Numericalsimulationsandvalidationsofsomesam-plesarethesubjectofongoingresearch.

Keywords: saturatedporousmedia;fluid-structurecoupling;principleofeffectivestress;dynamicmodelsofporosityandpermeability;consolidationequation

*

收稿日期: 2002-12-02

基金项目: 国家自然科学基金资助项目(10102020);宁波)中国科学院高新技术种子资金项目

: 1976~,男,,

420

水 动 力 学 研 究 与 进 展 2003年第4期

三维固结理论开展的,只是假设的多孔介质的应力应变本构关系有所不同(比如有的假设多孔介质为

1 引言

经典渗流力学一般假定流体流动的多孔介质(比如岩石、土壤等)是完全刚性的,即在孔隙流体压力变化过程中,固体骨架不产生任何弹性或者塑性变形,这时可将渗流视为非耦合问题来研究。这种简化可以得到问题的近似解,对于初期或者说以往的学科发展、工程实践都产生了巨大的和积极的作用,但是它存在着很多的缺陷。实际的多孔介质,不论是天然地质材料还是人造多孔固体,大多为可变形体,在实际的渗流过程中,由于孔隙流体压力的变化,一方面要引起多孔介质骨架有效应力变化,由此导致储层特性比如渗透率、孔隙度等的变化;另一方面,这些变化又反过来影响孔隙流体的流动和压力的分布。因此,在许多情况下,必需考虑孔隙流体在多孔介质中的流动规律及其对多孔介质本身的变形或者强度造成的影响,即考虑多孔介质内应力场与渗流场之间的相互耦合作用。

近年来,流固耦合渗流问题越来越受到人们的重视,这方面的研究涉及到很多工程领域,诸如石油、天然气开发过程中的流固耦合问题,地下水抽放和油气开采引起的地表沉降的流固耦合,软土地基固结沉降中的流固耦合问题,煤层瓦斯的耦合渗流和突出,地热开发理论中的流固耦合,核废料处理,水库诱发地震的流固耦合,岩波和坝基的稳定性问题等等。此项研究已经成为目前科学研究的热门课题。国内这方面的研究还处于起部阶段,有待进一步加强[2]。

本文利用基于多孔介质的有效应力原理,在考虑流体的可压缩性和多孔介质变形特性的前提下,利用平衡条件建立起完整的多孔介质流固耦合的渗流理论,并对模型进行简化处理和对比论证,这对于进一步更加准确地研究可变形多孔介质中变形和流体的流动特性有着重要的意义。

[1]

弹塑性,有的为粘弹性等等)或者孔隙流体假设为多相流体或单相流体的差别(有的称为非饱和和多孔介质或者饱和多孔介质)。

综上所述,目前所应用的流固耦合渗流模型几乎都是建立在Terzaghi有效应力原理之上的。而Terzaghi有效应力原理有其不足之处,本文作者从渗流力学的观点出发,在前人研究的基础上,推导出了基于多孔介质的有效应力原理。该有效应力原理包含了多孔介质的结构参数)孔隙度U,从而代替了以前诸多有效应力公式中广为使用的经验参数(各个模型的经验参数范围各异,取决于其适用对象和使用条件,比如比较常用的Biot常数[5,6,7]等);而且该有效应力公式是解析的,理论上是严格的,它可以适用于各种不同的多孔介质,而无须区分不同的适用对象和特定条件。

基于多孔介质的有效应力原理为:

Rij=(1-U)RspDij+Uij

s

[9,10]

(1)

其中Rij为固体颗粒之间的应力,p为孔隙流体的压力。令(1-U)Rij=Rij为固体骨架的有效应力

(Terzaghi有效应力)。则Terzaghi有效应力为

cRpDij=Rij-Uijs

c

(2)

该原理恰好体现了多孔介质中的流固耦合效

应,或者说是骨架变形和流体流动的耦合。其中,p为孔隙流体平均压力,对于饱和多孔介质而言,就是单相流体(比如对于饱和土而言,就是孔隙水)的压

c力;Rij为多孔介质固体骨架有效应力;U为多孔介质的孔隙度。多孔介质的各种机制的变形(骨架或者颗粒的弹、塑性变形,及蠕变等等),其结果都是导致孔隙度的变化,因此孔隙度可作为表征多孔介质结构状态(比如变形等)的重要特性参数。下一节:孔隙度U和渗透率K的动态模型即详细论述了孔隙度变化与多孔介质体积变形的关系。

2 有效应力原理

Terzaghi在研究土力学时,提出了其著名的有效应力原理[3],并建立了一维固结模型;之后Biot在此基础上取得了一些开创性的研究成

果,建立了较为完善的三维固结理论。从此以后,流[6,7,8,11][4,5]

3 孔隙度U和渗透率K的动态模型

对于流固耦合渗流,如引言所及,孔隙度U和渗透率K等物性参数当然是动态变化的,因此在建,

李培超等:饱和多孔介质流固耦合渗流的数学模型

421

诸多流固耦合渗流文献[6,7,8]所采用或者建立的模型中都忽略了这个重要的事实,即没有考虑孔隙度和渗透率的动态变化,这可能因为是受到经典渗流力学的影响疏忽所致。在经典渗流力学中,认为固体骨架不产生任何弹性或者塑性变形,自然不必考虑孔隙度和渗流率的动态变化。

文献[11]中提出,实现流固耦合渗流数值模拟的关键问题之一是如何建立流固耦合作用下的物性参数动态模型;并给出了用于计算孔隙度和渗透率的动态模型,但是该模型仍有些不尽合理的地方。

我们在文献[11]的基础上进行如下分析推导:

多孔介质固体骨架(Skeleton)体积用Vs表示,其变化我们用$Vs表示;多孔介质外观体积(BulkVolume)用Vb表示,其变化用$Vb表示;多孔介质孔隙体积(PoreVolume)用Vp表示,其变化用$Vp表示。

根据孔隙度的定义,我们有

Vp0+$VpVs0+$Vs

U==1-=

Vb0+$VbVb0+$Vb

Vs0(1+$Vs/Vs0)

=

Vb0(1+$Vb/Vb0)(1-U0)

(1+$Vs/Vs0)

1+EV

U=

U0+Ep/KsV+(1-U0)$

1+EV

(7)

1-U=1-0(1+$Vs/Vs0)=

1+EV

$Vs/Vs0=-$p/Ks+BTs$(5)式代入(3)式,有

(5)

$Vs/Vs0=BTs$

综合以上两者因素,我们有

0(1-$p/Ks+BT)s$1+EV

(6)

这里,$p=p-p0,$T=T-T0。

假定地层为等温,那么有

1-

U0+EV

。

1+EV

分析此式,如果存在U=U0,那么只有两种情况:U0=1(可认为是纯流体)或者E除V=0(无体应变)。此之外,没有U=U0的出现。

在等温的情形下,文献[11]给出U=

(3)

而实际上这并不是完全符合物理实际的[10]。分析式(7),除了上述的U0=1外,如果有-$p/Ks=EV,那么也可以有U=U0。

这种情况实际上就是同步变形的情形。也就是

(4)

说固体骨架的体积应变等于多孔介质总的体积应变,也等于孔隙体积的应变,即$Vs/Vs=$Vb/Vb=$Vp/Vp。图形可见图1。

1-

如果假定$Vs=0代入上式,可得到

U(1-U0+EV0)

U=1-=

1+E1+EVV

这与文献[11]的结论是一致的,然而实际上,即使不考虑温度场的影响,固体颗粒本身的体积也是

变化的,即肯定有本体变形[10],通常固体颗粒可认为是弹性体,所以本体变形通常是弹性变形,当然也可以是弹塑性变形。

如果认为固体颗粒是弹性变形,下边考虑孔隙流体压力对颗粒体积变形的影响。

因孔隙流体压力变化而引起的固体颗粒体积变化为

$Vs/Vs0=-$p/Ks

这种情形在自然界中显然是存在的,而文献

[11]并不能解释此现象。

再者,对于EV=0,代入其表达式,有U=U0+U0)p/s[图1 骨架的本体变形

422

水 动 力 学 研 究 与 进 展 2003年第4期

有-$p/Ks=0,即认为固体颗粒为刚体,没有任何压缩变形,那么才有U=U0。

在经典渗流力学中,通常认为固体骨架不变形,这一点可以接受。但是现在我们研究的是变形多孔介质中的流固耦合渗流问题,显然变形是研究的核心之一。所以如果忽略固体颗粒本身的变形是不妥当的。而且在工程实际中,固体颗粒总会或多或少有变形的。

根据以上推导和论述,本文提出的U动态模型是有效的而且是较为严格的。而文献[11]的动态模型存在一些缺陷,它不能解释一些物理现象,因此并不具有广谱性。

对于我们的公式,体应变EV可以是弹性,可以是塑性。通常状况下,它是弹塑性变形或者粘弹性变形等。而固体颗粒的体积变形,也可弹可塑,这里我们假定为弹性变形,而这也通常是符合物理实际的。

另外一种方法推导孔隙度动态模型。根据流体力学的连续性方程,可以得到固体骨架的连续性方程为(参考本文/5渗流场0部分的推导)

5Qs-Q=0(8)s考虑温度场效应后,上式可改写为:

(1-U0)(1-$p/Ks+Bs$T)

(13)

1+EV

U=1-

这与由孔隙度的定义出发推导出的孔隙度动态模型是完全一致的,从而再一次论证了本文孔隙度动态

模型的正确性。

同理,多孔介质绝对渗透率k并不是常数,在流固耦合渗流过程中,是受诸多因素制约并不断变化的,根据渗流力学Kozeny方程和文献[11],可以推导出k的方程为:

k0

#

1+EV

k=

E(Bs$T+$p/Ks)(1-U0)3V

[1+-](14)00

4 应力场方程

多孔介质有效应力应变本构关系方程为:

Rij=DijklEkl

c

Q)ý#V )s(1-Us+(1-U处理后得到

(15)

5EV(1-U)(+K)=s

对以上方程积分得到:

(1-U)=C#e

-(E+$p/K)

V

s

为描述简单起见,假定介质固体骨架为各向同

(9)

性弹性体(也可以是弹塑性本构关系,或者其它类型本构关系),Wx,Wy,Wz为固体骨架三个方向的位移,则有

cREEij=KVDij+2Lij

(10)

(16)

如果EV=0,Ks->+],那么USU0成立,代入上式,确定出C=1-U0。

至此,得到

U=U(EV,p)=1-(1-U0)#e将上式的指数项进行近似展开,得到:

(1-Up/Ks)0)(1-$

1+EV

(-$p/K-E)

s

V

其中EV为骨架体积应变,E =V=ý#W

5Wi5Wx5Wy5Wz

=++(17)i(11)

几何方程:

Eij=

(Wj,i+Wi,j)2

(18)

U=1-

(12)

应力平衡方程:

李培超等:饱和多孔介质流固耦合渗流的数学模型

423

Rij,j+Fi=0

联立式(2),式(19)有

c

RpDij,j+(Uij),j+Fi=0

(19)

; Vr为流体相对于骨架颗粒的速度,其表达式为 Vr=#

j

VjD,其中,下标j代表流体相, VjD为j相流体Dar-粒运动的绝对速度,根据定义有 Vs=cy速度,根据Darcy定律,我们有

KKrj

ý(pj-QjgD)j

(20)

将式(16)代入上式,我们便可以得到以基本未知量Wx,Wy,Wz,U,p为变量的平衡方程:

5EV+Lý2Wx+U=05EV

+Lý2Wy+U=0VjD=-

(21)

(28)

(K+L)

其中K为多孔介质的绝对渗透率张量,Krj为j相流体的相对渗透率。

对于本文,因为假定为单相流体(饱合多孔介

质),所以有Sj=1,进而有

Vr=-

ý(p-QfgD)(29)

(K+L)(22)

5EV2(K+L)+LýWz+U+fz=0(23)其中,fz=[(1-U)Qs+UQw]#g。对于各向同性弹性体,有

K=,G=L=,

(1+M)(1-2M)2(1+M)

则K+L=,代入式(21)~式(23)有

(1-2M)

5EV+Gý2Wx+U=0(24)

(1-2M)5EV2+GýWy+U(1-2M)=0(25)5EV2+GýWz+U+fz=0(26)

(1-2M)式(24)~式(26)即应力场方程,它对应于Biot[4]的三维固结方程(即其论文中的式(4.1)。

多孔介质骨架的连续性方程为

5[Qs(1-U)]

=0(30)

5t

ý#[Q) Vs]+s(1-U

孔隙流体(不考虑源汇项)的连续性方程为:

5[QfU]

=0ý#[QVf]+ fU

化简以上两式,得到:

(31)

Q)ý# Vs+(1-U)s(1-U

5Qs-Q=0s5t5t

(32)5Qf=0+Qf(33)

QVr+QVs+UfUý# fUý#

5 渗流场方程

由于渗流发生在可变形的多孔介质中,因而不但有流体具有一定的渗流速度,而且骨架颗粒也有一定的运动速度,所以流体质点的速度为

Vf= Vr+ Vs

(27)

将上两式两边分别除以Qs、Qf,并相加得到:Uý# Vr+ý#V s+

sf5Q5Q

+=0

Qs5tQf5t

(34)

通常而言,流体在等温条件下的状态方程可以

采用下式表达:

fs

424

[(p-Qf=Q0e

p0)/Kf]

水 动 力 学 研 究 与 进 展 2003年第4期

(35)

p|t=0=pi

(39)

其中Kf为孔隙流体的体积弹性压缩模量。所以有

5Qf=Kff同理,对于骨架固体颗粒,我们也可以有

p|边界=p1

s5Q=Kss其中pi是原始地层孔隙压力,即开采前孔隙流

体的压力分布。

(36)

(2)边界条件

渗流场边界通常有两种类型。定压边界:

(40a)

(37)

定流量(产量)边界条件:

(ýp-Qn|边界=qwgýD)# Ks为多孔介质骨架固体颗粒的体积弹性压缩模量。

V5E

)==,将式(36),式(37),式(28)代入式(34)得到多孔介质流固耦合渗流场方程(孔隙流体压力满足的微分

(40b)

而且ý# Vs=ý#(

方程):

5EV-ý#((ýp-QfgýD))++

(K+K)=0

sf

n为边界的法向量,q为流量或者产量。

如果是封闭边界(不渗透边界),则有(ýp

-QwgýD)=0。

应力模型边界条件有应力边界条件和位移边界条件。

位移边界条件:

W |边界=W 1为边界上的位移)(41)1 (W

(38)

应力边界条件:

上式中多孔介质绝对渗透率k并不是常数,具体用式(14)描述。

式(6),式(24)~式(26)和式(38)就构成了描

述单相流体的流固耦合渗流的数学模型。该模型共有5个方程,求解变量为Wx,Wy,Wz,U,p,总计5个,可见控制方程组是封闭的。

c

Rnj|边界=Ti (Ti为应力边界上的面力)ij#

(42)

式(6),式(24)~式(26)和式(38)以及定解条件式(39)~式(42)便构成了完整的单相流体流固耦合渗流的数学模型。对于此数学模型,我们可以采用有限差分、有限元等数值方法进行离散求解,从而求解出渗流场(p)的分布,孔隙度U以及位移场(Wi)的分布,进而更好地应用于油藏工程,岩土工程,地下水处理等领域,为生产实践提供一定的理论基础指导以及相应的工程资料。

该数学模型的数值求解和工程实例计算正在进行之中,可参考后继的相关研究和报道。

6 定解条件

对于上述控制方程组,应该补充以适当的边界条件和初始条件,才能构成定解问题。由于该模型包含渗流场方程和应力场方程,所以必需提供各自对应的边界条件。对于应力场方程,根据弹塑性力学基础知识,应该提供应力边界和位移边界;而对于渗流场方程,应该提供渗流边界。

(1)初始条件

通常指初始时刻或者从某一时刻起多孔介质孔7 模型的简化和定性分析

经典的Biot饱和土体一维竖向固结方程[4]为:

22=5zCV其中a=

CcV是竖向固结常数,

CcV

=K/(La),

(43)

在固结沉降过程中,孔隙度U和渗透率k是不断变化的,根据式(48),固结系数CV显然并非常量,这充分体现了流固耦合效应)变形与渗流的耦

合。它考虑了土体的重要物性参数孔隙度U的影响,并且孔隙流体以及固体颗粒的变形效应也被考虑进来。

在Biot固结理论中,假定孔隙流体和固体颗粒都是不可压缩的,而在我们的公式中,Ct显然包含孔隙流体和固体颗粒的压缩(分别以Kf,Ks表示)。如果不考虑孔隙流体和固体颗粒的压缩性,那么有UCt,进而式(48)变成:

=--Ua。

2(1-M)G

上述方程辅助以适当的初始和边界条件就构成了定解的一维固结沉降问题。

对于我们的模型而言,将方程式(24)~式(26)分别对x,y,z求导,相加得到:

22(1-2M+G)ýEV+Uýp=0将上式对时间t求导,并有a=,

2(1-M)G则

5EV=-Ua将式(44)代入式(38),有

-ý#((ýp-QfgýD))++-Ua)KsKf=0

CV=(49)

可见我们的固结系数比Biot固结系数多出因

(44)

,究其根源,在于Biot固结理论的根基是c

(广义的)Terzaghi有效应力原理:Rij=Rij-子-ApDij,其中A为Biot常数,A=1-Ks,K==3(1-2M)3(1-2M)。

而我们提出的基于多孔介质的有效应力原理

c

为:RpDij=Rij-Uij。因为基本原理的不同,从而导

((45)

致Biot三维固结方程[4](即其文章的式(4.1))和我们的控制方程式(24)~式(26)的差别,即前者中的-A被后者方程中的U代替,进而导致了固结系数的差别。

对于饱和土体,Biot认为A=1,A取为1实际上对应U趋近于1,这一点是有其物理含义的。根据孔隙度的定义,对于孔隙度U很高的情形,自然意味着孔隙体积较大,同时骨架体积较小或者说固体颗粒较少,那么多孔介质体积模量K与固体颗粒的体积模量Ks的比值变会较小。如果Uy1,那么,则Ay1。从而基Ksc

于多孔介质的有效应力原理RpDij=Rij-Uij就退化就有K/Ksy0,再由A=1-c为Rij=Rij-pDij,即Terzaghi有效应力原理的原

则上式相应的一维竖向固结控制方程为:2-+-Ua)=0(46)2(L5zKsKf5t令UCt=+,Ct是多孔介质总体压

KsKf

缩系数,方程变成

22=CV5z

这里,

t-Ua)

(47)

始形式;同时,CV也就变成了CcV,即此时我们的固

CV=

(48)

结方程退化成了经典的Biot固结方程。

所以可以认为Biot固结理论是我们提出的模型的特殊情形或者说是简化,同时这也一定程度上验证了我们所提出模型的正确性和可行性。

NewYork:Wiley,1943.

[4] BIOTMA.Generaltheoryofthreedimensionalconsol-i

dation[J].J.Appl.Phys1,1941,12:155-164.[5] BIOTMA.WillisDG.Theelasticcoefficientsofthe

theoryofconsolidation[J].ASMEJ.ApplMech.,1957,24:594-601.

[6] ZIENKIEWICZOCandSHIOMIT.Dynamicbehavior

ofsaturatedporousmedia:thegeneralizedBiotformula-tionanditsnumericalsolution[J].Int.J.Num.andAnaly.Meth.inGeomech,1984,8:71-96.

[7] CHENHY.Coupledfluidflowandgeomechanicsin

reservoirstudy-1.Theoryandgoverningequations[R].SPE30752,1995.

[8] 董平川,徐小荷.储层流固耦合的数学模型及其有限元

方程[J].石油学报,1998,19(1):64-70.

8 结论

本文经过数学推导,得到了饱合多孔介质流固耦合渗流数学模型,并对其进行简化分析和对比验证。该模型首先引入了基于多孔介质的有效应力原

理以代替通常使用的Terzaghi有效应力原理;孔隙度、渗透率动态模型也被引入控制方程中,这些都充分体现了流固耦合的物理特性,经离散后可用于流固耦合问题的数值计算。

参 考 文 献:

[1] 郭尚平,刘慈群,阎庆来.渗流力学的近况与展望[J].

力学与实践,1981,3(3):2-6.

[2] 董平川,徐小荷,何顺利.流固耦合问题及研究进展

[J].地质力学学报,1999,5(1):17-26.[3]

TERZAGHIK,TheorericalSoilMechanics[M].

[9] 徐献芝,李培超,李传亮.多孔介质有效应力原理研究

[J].力学与实践,2001,23(4):42-45.

[10]李传亮,孔祥言,徐献芝,李培超.多孔介质的双重有效

应力[J].自然杂志,1999,21(5):288-297.

[11]冉启全,李士伦.流固耦合油藏数值模拟中物性参数动

态模型研究[J].石油勘探与开发,1997,24(3):61-65.

本文来源：https://www.bwwdw.com/article/uez1.html