一种多孔粘弹性介质热渗流数值模拟方法与流程

1.本发明涉及地球动力学数值模拟技术领域，特别涉及一种多孔粘弹性介质热渗流数值模拟方法。


背景技术：

2.现有的多孔介质热渗流数值模拟程序模拟的介质模型多为二维平面；
3.且现有的多孔介质热渗流数值模拟程序耦合的的物理场不够全面，多为流
‑
固耦合和热
‑
固耦合；
4.现有的多孔介质热渗流数值模拟程序涉及的物性参数为常数，不符合实际多孔介质中物性参数随温度、压力等发生变化的规律；
5.现有的多孔介质热渗流数值模拟应用受限，只能数值模拟特定的问题；
6.现有的多孔介质热渗流数值模拟求解位移、应力场时，未考虑长时间尺度下模型介质表现为粘性，不符合实际情况；
7.现有的多孔介质热渗流数值模拟未考虑达西速度对热液在多孔介质中运移的影响；
8.现有的多孔介质热渗流数值模拟程序不具有扩展性，不能引入其它物理场，不能改进算法等。


技术实现要素：

9.为了至少解决或部分解决上述问题，提供一种多孔粘弹性介质热渗流数值模拟方法。
10.为了达到上述目的，本发明提供了如下的技术方案：
11.本发明一种多孔粘弹性介质热渗流数值模拟方法，包括以下步骤：
12.s1：根据初始条件、边界条件、温度、温度变化速率、体应力、有效应力，通过孔隙渗流扩散模型得到孔隙水压，并计算出孔隙水压变化速率，由达西定律计算出达西速度，将孔隙水压、达西速度值传递至孔隙介质热传递模型，将孔隙水压、孔隙水压变化速率值传递至粘弹性体本构模型；
13.s2：根据初始条件、边界条件、孔隙水压、达西速度，通过孔隙介质热传递模型得到温度值，并计算出温度变化速率，将温度、温度变化速率传递至粘弹性体本构模型和孔隙渗流扩散模型；
14.s3：根据初始条件、边界条件、温度、温度变化速率、孔隙水压、孔隙水压变化速率，求解出粘弹性体位移、应力数值解，并计算出有有效应力，将体应力、有效应力值传递至孔隙渗流扩散模型，将有效应力值传递至孔隙介质热传递模型；
15.s4：重复步骤s1
‑
s3，直至计算结束。
16.作为本发明的一种优选技术方案，所述孔隙渗流扩散模型的计算公式为
[0017][0018]
其中，σ
kk
为体应力，b为skempton系数，p为孔隙水压，t为温度，k为渗透系数，μ为粘滞系数，g为剪切模量，υ为泊松比，υ
u
为不排水泊松比，α为热膨胀系数。
[0019]
作为本发明的一种优选技术方案，所述孔隙度经验公式为
[0020][0021]
其中，为孔隙度，为初始孔隙度，α为实验系数，p
e
为有效应力。
[0022]
作为本发明的一种优选技术方案，所述孔隙介质热传递模型的计算公式为
[0023][0024]
其中，为孔隙度，ρ
f
为热液密度，ρ
s
为固体密度，c
f
为热液比热容，c
s
为固体比热容，t为温度，v为达西速度，λ
f
为热液热传导率，λ
s
为固体热传导率。
[0025]
作为本发明的一种优选技术方案，所述渗透率经验公式为
[0026][0027]
其中，k为渗透率，k0为初始渗透率，b为实验系数，p
e
为有效应力。
[0028]
作为本发明的一种优选技术方案，所述粘弹性本构模型的计算公式为
[0029][0030]
其中：e、g、v、λ、μ分别为包含了孔隙、裂隙在内的总体材料的宏观等效弹性模量、剪切模量、泊松比、拉梅常数，υ
u
为不排水泊松比，σ
ij
为应力矢量，ε
ij
为应变矢量，b为skempton系数，p为孔隙水压，为孔隙水压变化速率，t为温度，为温度变化速率，α为热膨胀系数，η为粘滞系数；
[0031][0032]
平衡方程：
[0033]
σ
ij，j
f
i
＝0
[0034]
几何方程：
[0035]
[0036]
作为本发明的一种优选技术方案，所述热液粘度经验公式为
[0037][0038]
其中，μ为粘度，t为温度。
[0039]
与现有技术相比，本发明的有益效果如下：
[0040]
本发明可以模拟二维和三维模型，将二维扩展到三维，具有更高的应用价值，可供科学与工程计算之用；
[0041]
本发明为流
‑
热
‑
固耦合的有限元程序，既符合实际物理现象，又可以定量的研究各物理场间的相互影响；
[0042]
本发明引入的物性参数为温度、有效应力的经验函数，符合实际物理现象。
[0043]
本发明应用范围较广，所模拟的多孔介质可以是天然的(如植物、石头、土壤等)、人造的(如砾石、混凝土、过滤器、凝胶等)、人体骨骼等，只需对物性参数、初始条件、边界条件稍作修改，就可以模拟其它多孔介质模型热渗流现象；
[0044]
本发明引入多孔介质粘弹性本构关系，通过改变粘度，可同时模拟脆性、韧性多孔介质模型；
[0045]
本发明引入达西速度和热对流，能够模拟热液受达西速度的影响下在多孔介质中的运移情况；
[0046]
本发明既可以在原有程序的基础上增加其他物理场，又可以改进算法，具有高度扩展的特性。
附图说明
[0047]
附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：
[0048]
图1是本发明的模型之间的数据传递图；
[0049]
图2是本发明验证孔隙渗流扩散以及孔隙介质热传递程序的坐标系图；
[0050]
图3是对比图2模型中心点(x＝50，y＝50)处解析解与数值解的结果图；
[0051]
图4是本发明验证maxwell体粘弹性本构关系部分的程序的坐标系图；
[0052]
图5是对比图4模型中心点(x＝15，y＝300)处解析解与数值解的结果图。
具体实施方式
[0053]
以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。
[0054]
此外，如果已知技术的详细描述对于示出本发明的特征是不必要的，则将其省略。
[0055]
实施例1
[0056]
如图1所示，本发明提供一种多孔粘弹性介质热渗流数值模拟方法，包括以下步骤：
[0057]
s1：根据初始条件、边界条件、温度、温度变化速率、体应力、有效应力，通过孔隙渗流扩散模型得到孔隙水压，并计算出孔隙水压变化速率，由达西定律计算出达西速度，将孔隙水压、达西速度值传递至孔隙介质热传递模型，将孔隙水压、孔隙水压变化速率值传递至
粘弹性体本构模型；
[0058]
s2：根据初始条件、边界条件、孔隙水压、达西速度，通过孔隙介质热传递模型得到温度值，并计算出温度变化速率，将温度、温度变化速率传递至粘弹性体本构模型和孔隙渗流扩散模型；
[0059]
s3：根据初始条件、边界条件、温度、温度变化速率、孔隙水压、孔隙水压变化速率，求解出粘弹性体位移、应力数值解，并计算出有有效应力，将体应力、有效应力值传递至孔隙渗流扩散模型，将有效应力值传递至孔隙介质热传递模型；
[0060]
s4：重复步骤s1
‑
s3，直至计算结束。
[0061]
该发明方法是依托有限元程序自动生成系统(finite element program generator，fepg)，将多孔介质热渗流微分模型以及算法文件编译成fortran源码，供有限元软件开发和科学与工程计算之用。
[0062]
多孔介质热渗流数值模拟涉及到三个微分模型：孔隙渗流扩散模型、孔隙介质热传递模型、(粘)弹性体本构模型。下面一一介绍各个微分：
[0063]
多孔介质热渗流数值模拟微分模型
[0064]
(一)孔隙渗流扩散模型
[0065][0066]
其中，σ
kk
为体应力，b为skempton系数，p为孔隙水压，t为温度，k为渗透系数，μ为粘滞系数，g为剪切模量，υ为泊松比，υ
u
为不排水泊松比，α为热膨胀系数。
[0067]
孔隙渗流扩散方程以固结理论为基础，引入温度的影响，描述了孔隙水压在多孔介质中的传递。
[0068]
流体在多孔介质中流动时,流体与多孔介质之间不可避免地要发生交互作用,即流—固耦合作用。多孔介质骨架在流体载荷(孔隙水压)的作用下会发生变形和运动。变形和运动又会反过来影响其中流体的流动,从而改变流体载荷(孔隙水压)的分布和大小。对于地下含水层,流体运动过程中不仅会发生流—固耦合现象,同时还受到外部载荷、温度的综合影响。
[0069]
(二)孔隙介质热传递模型
[0070][0071]
其中，为孔隙度，ρ
f
为热液密度，ρ
s
为固体密度，c
f
为热液比热容，c
s
为固体比热容，t为温度，v为达西速度，λ
f
为热液热传导率，λ
s
为固体热传导率；
[0072]
对于瞬态问题流体的热传递需要考虑固相和液相各自的热传递及固液相之间的热交换。但是，如果流体和固体之间局部传热的响应时间比相位平均温度变化的特征时间小几个数量级，那么流体和固相温度可以假定为相同，即每个点局部热平衡(the local thermal equilibrium,缩写为lte)成立。在lte的假设下，仅需要单个能量方程来描述多孔介质内的热传递过程。
[0073]
假定流体和固体介质具有相同的温度t,在整个介质里都会发生传导热的渗透，并且方程右边的热传导率的适当值必须是整个充满流体的多孔介质和固体介质上的体平均值λ
m
。
[0074]
由于介质的基本成分是由此流体导热更好的固体介质组成的，因此假定λ
m
为固体介质的热导率一般都是一种很好的近似。热能被贮存于充满流体的孔隙和固体介质中。因此，方程中左边的热惯量项也是体平均值。由于只有流体传输热量，方程左端的对流项使用流体密度ρ
f
和流体比热c
f
。
[0075]
(三)粘弹性本构模型
[0076][0077]
其中：e、g、v、λ、μ分别为包含了孔隙、裂隙在内的总体材料的宏观等效弹性模量、剪切模量、泊松比、拉梅常数，υ
u
为不排水泊松比，σ
ij
为应力矢量，ε
ij
为应变矢量，b为skempton系数，p为孔隙水压，为孔隙水压变化速率，t为温度，为温度变化速率，α为热膨胀系数，η为粘滞系数；
[0078][0079]
平衡方程：
[0080]
σ
ij，j
f
i
＝0
[0081]
几何方程：
[0082][0083]
粘弹性介质选用maxwell体模型，在短时间尺度内表现为弹性体，在长的时间尺度则表现为流体，通过调整粘滞系数可以模拟脆性介质和柔性介质的变形。由于maxwell体具有这种可以较自然处理脆
‑
韧性转变的优点，利用三维maxwell体模型，使用有限元方法，通过解粘弹性本构模型、平衡方程、几何方程计算研究区域内的位移场和应力场。
[0084]
(四)物性参数经验公式
[0085][0086]
其中，为孔隙度，为初始孔隙度，a为实验系数，p
e
为有效应力；
[0087]
渗透率经验公式为
[0088][0089]
其中，k为渗透率，k0为初始渗透率，b为实验系数，p
e
为有效应力；
[0090]
热液粘度经验公式为
[0091][0092]
其中，μ为粘度，t为温度。
[0093]
多孔介质热渗流数值模拟程序是将上述微分方程通过有限元自动生成系统编译生成的有限元程序。
[0094]
程序正确性验证
[0095]
(1)为了验证孔隙渗流扩散以及孔隙介质热传递程序的正确性，选取一个有解析解(u＝cost(x2 y2))的算例。模型及坐标系选取如图2所示，块体的高h＝100m，宽l＝100m，从t＝0时刻起，计算时间10s，计算结果如图3所示。
[0096][0097]
(2)为了验证maxwell体粘弹性本构关系部分的程序的正确性，设计了一个有解析解的算例。
[0098]
模型及坐标系选取如图4所示，块体的高h＝600m，宽l＝30m，从t＝0时刻起在其顶面施加p0＝105pa，弹性模量e＝1.0
×
10
10
pa，泊松比v＝0.25，粘性系数η＝1.0
×
10
10
pa
·
s，计算时间25s；计算结果如图5所示。
[0099]
最后应说明的是：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。