CN109900610A - 一种多孔材料的联合测试方法 - Google Patents

Abstract

本发明属于多孔材料表征技术领域，具体涉及一种岩石等多孔低渗材料的绝对渗透率、滑移因子和有效孔隙率的联合测试方法。本发明包括如下步骤：步骤1、搭建试验系统；步骤2、实现非稳态渗流过程；步骤3、实现稳态渗流过程；步骤4、构建非稳态气体渗流控制方程；步骤5、利用全隐式有限体积格式求解非稳态气体渗流控制方程；步骤6、基于稳定渗流结果计算绝对渗透率k_∞与滑移因子b的对应关系；步骤7、基于最小二乘原理计算滑移因子、绝对渗透率和有效孔隙率。本发明能够解决目前绝对渗透率、滑移因子和有效孔隙率测试效率低下、理论基础不协调的问题。

Description

一种多孔材料的联合测试方法

技术领域

本发明属于多孔材料表征技术领域，具体涉及一种岩石等多孔低渗材料的绝对渗透率、滑移因子和有效孔隙率的联合测试方法。

背景技术

在页岩气开采、地下储气库开发以及放射性废物处置库研发等能源工程领域，硐室围岩以及混凝土等工程材料的气相绝对渗透率、滑移因子和有效孔隙率是进行开采效率预测、密封性评估和安全性评价必需的核心参数。绝对渗透率和滑移因子反映多孔材料中气体渗流通过的能力，绝对渗透率仅与孔隙结构相关，滑移因子则受到气体类型和孔隙结构的控制。有效孔隙率是指多孔材料中作为气体渗流通道的孔隙所占材料总体积的比例。然而，上述天然以及人工材料通常为低渗材料，其渗透率和孔隙率极小，测试较为困难。建立高效可靠的绝对渗透率、滑移因子和有效孔隙率测量方法，是提高天然和人工多孔材料性能评价可靠性的要求，对上述能源工程具有重要意义。

目前，尚未出现低渗材料的绝对渗透率、滑移因子和有效孔隙率的联合测试方法，需要不同的物理试验过程分别测量。渗透率的测试方法通常有稳态法和脉冲法两类，绝对渗透率和滑移因子则需要对不同压力条件下渗透率测试结果进行拟合计算获得。由于需要进行多次测试，整个试验的效率较低。

有效孔隙率的测定主要有流体置换法和直接扫描法两类。流体置换法在测定有效孔隙率时，无法对连通孔隙和开口不连通孔隙进行区分。由于包含了不参与渗流的开口不连通孔隙，上述方法测试获得的“有效孔隙率”比真正参与渗流的有效孔隙率要大。而CT等扫描法主要是对材料进行二维截面采样，继而进行三维建模，获得虚拟的三维孔隙结构，计算获得的有效孔隙率的精准度较低，同时，代价巨大，测试效率较低。总之，目前常规方法测试获得的有效孔隙率并不能准确反映控制多孔材料渗透特性的孔隙率。

因此，在多孔材料性能评价同时使用上述来自于不同物理过程的参数，理论基础并不协调，进而可能影响上述能源工程安全评价的可靠性。

发明内容

本发明解决的技术问题：本发明提供一种多孔材料绝对渗透率、滑移因子和有效孔隙率的联合测试方法，解决目前绝对渗透率、滑移因子和有效孔隙率测试效率低下、理论基础不协调的问题。

本发明采用的技术方案：

一种多孔材料的联合测试方法，包括如下步骤：

步骤1、搭建试验系统；

本试验系统包括气压加载装置、围压加载装置和试验腔体；圆柱形试样放置于圆形试验腔体中部，上下端分别与出气孔和进气孔相连，出气孔直接通向大气，进气孔连接气压加载装置，并在连接管路上安装阀门；围压加载装置连接到试验腔体，给试样施加相当于最大进气压2～5倍的围压；

步骤2、实现非稳态渗流过程；

步骤3、实现稳态渗流过程；

步骤4、构建非稳态气体渗流控制方程；

步骤5、利用全隐式有限体积格式求解非稳态气体渗流控制方程；

步骤6、基于稳态渗流结果计算绝对渗透率k_∞与滑移因子b的对应关系；

步骤7、基于最小二乘原理计算滑移因子、绝对渗透率和有效孔隙率；

所述步骤1中，气压加载装置需要能够实现固定体积下的气压监测记录以及固定压力下的体积变化监测记录。

所述步骤1中，围压加载器可采用全自动或手动液压加载装置。

所述步骤2的具体步骤为，关闭气压加载器和进气孔之间的阀门，在气压加载器中形成气压P₀，即P(x＝0,t＝0)＝P₀，形成的储气空间体积为V(t)≡V₀；

之后打开阀门，在压力差作用下气体由试样进气端渗流向出气端；同时，储气空间中的气体摩尔数由于渗流而逐渐减小，因此，储气空间中的气压逐渐下降；记录气压下降过程，形成P(0,t)v.s.t曲线；

当储气空间中的气压降到P₁时，结束非稳定渗流过程，将气压加载器调整为稳压加载模式。

所述P₁为小于初始压力P₀的非稳态渗流结束控制压力，取值为0.8～0.9P₀。

所述步骤3的具体步骤为：进气端气压P₁保持不变，记录稳压加载过程中的储气空间体积变化过程，形成V(t)v.s.t曲线；实时计算曲线斜率，当两次计算相对差异不大于ζ时，取当前斜率的绝对值作为P₁压力条件下的气体流量Q。

所述ζ取值为1％。

所述步骤4的具体步骤为：

多孔材料渗流质量守恒方程：

多孔材料渗流达西定律：

多孔材料渗流气体滑移方程：

其中，φ为多孔材料孔隙率；t为时间；P为多孔材料孔隙中的气体压力，通常为空间、时间以及温度的函数；v为等效流速，此流速并不是气体在孔隙中的真实流动速度，而是以试样整个截面面积作为流动截面的等效流速；Q为气体流量；A为试样垂直于渗流方向上的截面面积；k为试样有效气体渗透率，为空间和时间的函数；μ为气体粘性系数，可以看作恒定常数；k_∞为绝对渗透率，仅与孔隙结构相关；b为滑移因子，与气体类型和孔隙结构相关；

根据公式(1)、(2)、(3)，并定义U(x,t)＝P(x,t)+b，获得气体渗流控制方程：

所述步骤5的具体步骤为：

非稳态渗流过程的初始状态可以表述为：

非稳态渗流过程的出气端边界条件可以表述为：

U(H,t)＝P_atm+b (6)

非稳态渗流过程的进气端边界条件，可以通过物质守恒原理获得：

其中，P_atm为大气压，可以取为1×10⁵Pa。

步骤4中获得控制方程以及上述初始状态和边界条件方程，组成一个二维非线性抛物线型偏微分方程定解问题；对于此类非线性问题，采用有限体积法中的全隐式Picard-Newton迭代格式进行数值求解；Picard-Newton迭代格式具有良好的二阶收敛特性；为构建Picard-Newton迭代格式，对控制方程左边时间项进行向后欧拉时间离散；对控制方程右边的k_∞U/μ项进行Taylor级数展开，继而完成Newton线性化，获得如下线性迭代格式：

式中，j为当前时间步标号，可以取为1,2,…,n-1；n为总时间步数；j＝0对应于初始状态；(s)和(s+1)分别代表第s次和第s+1次迭代计算结果，s可以取为0,1,2,…；Δt为时间步长(s)；

采用全隐Picard空间离散格式，则上式最终可以变换为：

式中：i为当前空间网格点编号，可以取为1,2,…,m-1；i＝0为进气端；i＝m为出气端；为第i个空间网格点在j-1时间步上的最终解，以下表达式中未带迭代计算次数(s)的U值均为最终解。

为启动公式(9)的计算，需要定义所有网格点的初始压力值；

代表初始条件的方程可以离散为：

公式(9)形成含有m+1个未知量的m-1个方程，需要带入上下边界的两个边界条件方程，才能完成计算；代表边界条件的方程可以离散为：

其中

同时，需要假定第j时间步的初始值为第j-1时间步的最终解，即：

根据上述离散的控制方程(公式(9))、初始条件方程(公式(10)和(11))和边界条件方程(公式(12)和(13))，可以形成m+1个线性方程，继而可以求解此线性方程组m+1个未知量以此类推可以求解获得当两次迭代计算结果差值的最大值小于精度要求时，即认为获得第j时间步的最终解。其中，ξ为气压最小分辨率，取为10^-8P₀。

所述步骤6的具体步骤为：

对于稳定渗流，步骤4控制方程左侧的时间域偏微分项为0，控制方程具有解析解：

将上式以及步骤2获得的P₁压力条件下的流量Q带入步骤4所示的渗流达西定律，可以获得绝对渗透率k_∞与滑移因子b的换算关系：

其中，H为试样高度(m)。

所述步骤7的具体步骤为：

假定滑移因子为b、有效孔隙率为φ，并根据公式(16)计算绝对渗透率k_∞；利用步骤5中公式(9)、(10)、(11)、(12)和(13)形成线性方程组，计算获得的第j时间步的储气空间压力为采用最小二乘原理设定目标函数如下：

式中，为步骤2获得的P(0,t)v.s.t曲线在j时间步处的实测压力值；采用不需要直接计算梯度的修正Powell方向加速法进行计算上述目标函数的最小值，继而获得绝对渗透率k_∞、滑移因子b和有效孔隙率φ最终结果。

为求解上述目标函数最小值及其对应的绝对渗透率k_∞、滑移因子b和有效孔隙率φ，采用不需要直接计算梯度的修正Powell方向加速法进行计算；主要计算步骤为：

步骤7.1：

定义算法起始点v₀ ⁽¹⁾＝[log(b_r),log(φ_r)]、初始搜索方向s₁ ⁽⁰⁾＝[1,0]和s₂ ⁽⁰⁾＝[0,1]、收敛精度ε＝10^-8。其中，b_r和φ_r分别为设定的滑移因子和有效孔隙率的初始值，根据测试材料的经验值进行确定，对于花岗岩可以分别取为1.0×10^-19m²和0.005；

步骤7.2：

对于第j轮计算(j＝1,2,…)，沿方向s₁ ^(j)和s₂ ^(j)逐次进行一维搜索，即

步骤7.3：

分别进行以下计算

判断下述两个条件是否同时满足

f₃＜f₁ (27)

步骤7.4：

当上述两个条件同时满足时，进行以下线性搜索

令

如果停止计算，最终参数值为

否则，重复步骤7.2和7.3进行第j+1轮计算；

步骤7.5：

当步骤7.3中两个条件不能同时满足时，则取

如果停止计算，最终参数值为

否则，重复步骤7.2和7.3进行第j+1轮计算；

步骤7.6：

根据步骤7.4或7.5获得的滑移因子b，利用公式(16)计算绝对渗透率k_∞。

本发明的有益效果：

本发明提供的方法测试过程简单，可在同一测试过程中获得绝对渗透率、滑移因子和有效孔隙率，并能够解决常规测试方法测试效率低下、理论基础不协调的问题。

附图说明

图1为本发明使用的测试系统组成示意图。

图2为本发明非稳态渗流路径示意图。

图3为本发明稳态渗流路径示意图。

图4为本发明非稳态渗流及稳态渗流过程中的进气端气压演化过程示意图。

图1中：1-1为气压加载装置、1-2为可变储气空间、1-3为围压加载装置、1-4为阀门、1-5为试验腔体、1-6为出气孔、1-7为进气孔、1-8为进气端连通管、1-9为下压头、1-10为液压油、1-11为试样、1-12为橡胶套、1-13为上压头、1-14为出气端连通管。

具体实施方式

下面结合附图和具体实施例对本发明提供一种多孔材料绝对渗透率、滑移因子和有效孔隙率的联合测试方法作进一步详细说明。

本发明提供一种多孔材料绝对渗透率、滑移因子和有效孔隙率的联合测试方法，包括如下步骤：

步骤1、搭建试验系统；

本发明使用的试验系统如图1所示。本试验系统包括气压加载装置(1-1)、围压加载装置(1-3)和试验腔体(1-5)。圆柱形试样(1-11)放置于圆形试验腔体(1-5)中部，上下端分别安装有上压头(1-13)和下压头(1-9)，并与出气孔(1-6)和进气孔(1-7)相连。出气孔(1-6)直接通向大气。进气孔(1-7)连接气压加载装置(1-1)，并在连接管路上安装阀门(1-4)。气压加载装置(1-1)需要能够实现固定体积下的气压监测记录以及固定压力下的体积变化监测记录，例如，可以采用英国GDS仪器设备有限公司生产的气体加压装置。围压加载装置(1-3)连接到试验腔体(1-5)，给试样(1-11)施加相当于最大进气压2～5倍的围压。围压加载器1-3可采用全自动或手动液压加载装置。

步骤2、实现非稳态渗流过程；

本发明形成的1维非稳态渗流路径如图2所示。

关闭气压加载器1-1和进气孔1-7之间的阀门1-4，在气压加载器1-1中形成气压P₀，即P(x＝0,t＝0)＝P₀，形成的储气空间体积为V(t)≡V₀。

之后打开阀门1-4，在压力差作用下气体由试样1-11进气端渗流向出气端。同时，储气空间中的气体摩尔数由于渗流而逐渐减小，因此，储气空间中的气压逐渐下降。记录气压下降过程，形成P(0,t)v.s.t曲线，如图4中非稳态渗流曲线部分所示。

当储气空间中的气压降到P₁时，结束非稳定渗流过程，将气压加载器调整为稳压加载模式。其中，P₁为小于初始压力P₀的非稳态渗流结束控制压力，可以取为0.8～0.9P₀。

步骤3、实现稳态渗流过程

本发明形成的1维稳态渗流路径如图3所示。

稳态渗流过程中，进气端气压P₁保持不变，如图4中稳态渗流曲线部分所示。记录稳压加载过程中的储气空间体积变化过程，形成V(t)v.s.t曲线。实时计算曲线斜率，当两次计算相对差异不大于ζ时，取当前斜率的绝对值作为P₁压力条件下的气体流量Q。其中，ζ为远小于1的误差控制指标，取为1％。

步骤4、构建非稳态气体渗流控制方程

多孔材料渗流质量守恒方程：

多孔材料渗流达西定律：

多孔材料渗流气体滑移方程：

其中，φ为多孔材料孔隙率(-)；t为时间(s)；P为多孔材料孔隙中的气体压力(Pa)，通常为空间、时间以及温度的函数，对于图2所示的1维恒温渗流过程，仅为坐标x和时间t的函数，即P(x,t)；v为等效流速(m/s)，此流速并不是气体在孔隙中的真实流动速度，而是以试样整个截面面积作为流动截面的等效流速；Q为气体流量(m³/s)；A为试样垂直于渗流方向上的截面面积(m²)；k为试样有效气体渗透率(m²)，为空间和时间的函数；μ为气体粘性系数(Pa·s)，可以看作恒定常数；k_∞为绝对渗透率(m²)，仅与孔隙结构相关；b为滑移因子(Pa)，与气体类型和孔隙结构相关。

根据上述多孔材料渗流质量守恒方程、达西定律和气体滑移方程，并定义U(x,t)＝P(x,t)+b，可以获得气体渗流控制方程：

步骤5、利用全隐式有限体积格式求解非稳态气体渗流控制方程

非稳态渗流过程的初始状态可以表述为：

非稳态渗流过程的出气端边界条件可以表述为：

U(H,t)＝P_atm+b (6)

非稳态渗流过程的进气端边界条件，可以通过物质守恒原理获得：

其中，P_atm为大气压，可以取为1×10⁵Pa。

步骤4中获得控制方程以及上述初始状态和边界条件方程，组成一个二维非线性抛物线型偏微分方程定解问题。对于此类非线性问题，采用有限体积法中的全隐式Picard-Newton迭代格式进行数值求解。Picard-Newton迭代格式具有良好的二阶收敛特性。为构建Picard-Newton迭代格式，对控制方程左边时间项进行向后欧拉时间离散；对控制方程右边的k_∞U/μ项进行Taylor级数展开，继而完成Newton线性化，获得如下线性迭代格式：

式中，j为当前时间步标号，可以取为1,2,…,n-1；n为总时间步数；j＝0对应于初始状态；(s)和(s+1)分别代表第s次和第s+1次迭代计算结果，s可以取为0,1,2,…；Δt为时间步长(s)。

采用全隐Picard空间离散格式，则上式最终可以变换为：

式中：i为当前空间网格点编号，可以取为1,2,…,m-1；i＝0为进气端；i＝m为出气端；为第i个空间网格点在j-1时间步上的最终解，以下表达式中未带迭代计算次数(s)的U值均为最终解。

为启动公式(9)的计算，需要定义所有网格点的初始压力值。代表初始条件的方程可以离散为：

公式(9)形成含有m+1个未知量的m-1个方程，需要带入上下边界的两个边界条件方程，才能完成计算。代表边界条件的方程可以离散为：

其中

同时，需要假定第j时间步的初始值为第j-1时间步的最终解，即：

根据上述离散的控制方程(公式(9))、初始条件方程(公式(10)和(11))和边界条件方程(公式(12)和(13))，可以形成m+1个线性方程，继而可以求解此线性方程组m+1个未知量以此类推可以求解获得当两次迭代计算结果差值的最大值小于精度要求时，即认为获得第j时间步的最终解。其中，ξ为气压最小分辨率，取为10^-8P₀。

步骤6、基于稳态渗流结果计算绝对渗透率k_∞与滑移因子b的对应关系

对于稳态渗流，步骤4控制方程左侧的时间域偏微分项为0，控制方程具有解析解：

将上式以及步骤2获得的P₁压力条件下的流量Q带入步骤4所示的渗流达西定律，可以获得绝对渗透率k_∞与滑移因子b的换算关系：

其中，H为试样高度(m)。

步骤7、基于最小二乘原理计算滑移因子、绝对渗透率和有效孔隙率

假定滑移因子为b、有效孔隙率为φ，并根据公式(16)计算绝对渗透率k_∞。利用步骤5中公式(9)、(10)、(11)、(12)和(13)形成线性方程组，计算获得的第j时间步的储气空间压力为采用最小二乘原理设定目标函数如下：

式中，为步骤2获得的P(0,t)v.s.t曲线在j时间步处的实测压力值。采用不需要直接计算梯度的修正Powell方向加速法进行计算上述目标函数的最小值，继而获得绝对渗透率k_∞、滑移因子b和有效孔隙率φ最终结果。主要计算步骤为：

步骤7.1：

定义算法起始点v₀ ⁽¹⁾＝[log(b_r),log(φ_r)]、初始搜索方向s₁ ⁽⁰⁾＝[1,0]和s₂ ⁽⁰⁾＝[0,1]、收敛精度ε＝10^-8。其中，b_r和φ_r分别为设定的滑移因子和有效孔隙率的初始值，根据测试材料的经验值进行确定，对于花岗岩可以分别取为1.0×10^-19m²和0.005。

步骤7.2：

对于第j轮计算(j＝1,2,…)，沿方向s₁ ^(j)和s₂ ^(j)逐次进行一维搜索，即

步骤7.3：

分别进行以下计算

判断下述两个条件是否同时满足

f₃＜f₁ (27)

步骤7.4：

当上述两个条件同时满足时，进行以下线性搜索

令

如果停止计算，最终参数值为

否则，重复步骤7.2和7.3进行第j+1轮计算。

步骤7.5：

当步骤7.3中两个条件不能同时满足时，则取

如果停止计算，最终参数值为

否则，重复步骤7.2和7.3进行第j+1轮计算。

步骤7.6：

根据步骤7.4或7.5获得的滑移因子b，利用公式(16)计算绝对渗透率k_∞。

Claims (13)

1.一种多孔材料的联合测试方法，其特征在于：包括如下步骤：

步骤(1)、搭建试验系统；

本试验系统包括气压加载装置、围压加载装置和试验腔体；圆柱形试样放置于圆形试验腔体中部，上下端分别与出气孔和进气孔相连，出气孔直接通向大气，进气孔连接气压加载装置，并在连接管路上安装阀门；围压加载装置连接到试验腔体，给试样施加相当于最大进气压2～5倍的围压；

步骤(2)、实现非稳态渗流过程；

步骤(3)、实现稳态渗流过程；

步骤(4)、构建非稳态气体渗流控制方程；

步骤(5)、利用全隐式有限体积格式求解非稳态气体渗流控制方程；

步骤(6)、基于稳态渗流结果计算绝对渗透率k_∞与滑移因子b的对应关系；

步骤(7)、基于最小二乘原理计算滑移因子、绝对渗透率和有效孔隙率。

2.根据权利要求1所述的一种多孔材料的联合测试方法，其特征在于：所述步骤(1)中，气压加载装置需要能够实现固定体积下的气压监测记录以及固定压力下的体积变化监测记录。

3.根据权利要求1所述的一种多孔材料的联合测试方法，其特征在于：所述步骤(1)中，围压加载器可采用全自动或手动液压加载装置。

4.根据权利要求3所述的一种多孔材料的联合测试方法，其特征在于：所述步骤(2)的具体步骤为，关闭气压加载器和进气孔之间的阀门，在气压加载器中形成气压P₀，即P(x＝0,t＝0)＝P₀，形成的储气空间体积为V(t)≡V₀；

之后打开阀门，在压力差作用下气体由试样进气端渗流向出气端；同时，储气空间中的气体摩尔数由于渗流而逐渐减小，因此，储气空间中的气压逐渐下降；记录气压下降过程，形成P(0,t)v.s.t曲线；

当储气空间中的气压降到P₁时，结束非稳定渗流过程，将气压加载器调整为稳压加载模式。

5.根据权利要求4所述的一种多孔材料的联合测试方法，其特征在于：所述P₁为小于初始压力P₀的非稳态渗流结束控制压力，取值为0.8～0.9P₀。

6.根据权利要求5所述的一种多孔材料的联合测试方法，其特征在于：所述步骤(3)的具体步骤为：进气端气压P₁保持不变，记录稳压加载过程中的储气空间体积变化过程，形成V(t)v.s.t曲线；实时计算曲线斜率，当两次计算相对差异不大于ζ时，取当前斜率的绝对值作为P₁压力条件下的气体流量Q。

7.根据权利要求6所述的一种多孔材料的联合测试方法，其特征在于：所述ζ取值为1％。

8.根据权利要求7所述的一种多孔材料的联合测试方法，其特征在于：所述步骤(4)的具体步骤为：

多孔材料渗流质量守恒方程：

多孔材料渗流达西定律：

多孔材料渗流气体滑移方程：

其中，φ为多孔材料孔隙率；t为时间；P为多孔材料孔隙中的气体压力，通常为空间、时间以及温度的函数；v为等效流速，此流速并不是气体在孔隙中的真实流动速度，而是以试样整个截面面积作为流动截面的等效流速；Q为气体流量；A为试样垂直于渗流方向上的截面面积；k为试样有效气体渗透率，为空间和时间的函数；μ为气体粘性系数，可以看作恒定常数；k_∞为绝对渗透率，仅与孔隙结构相关；b为滑移因子，与气体类型和孔隙结构相关；

根据公式(1)、(2)、(3)，并定义U(x,t)＝P(x,t)+b，获得气体渗流控制方程：

9.根据权利要求8所述的一种多孔材料的联合测试方法，其特征在于：所述步骤(5)的具体步骤为：

非稳态渗流过程的初始状态可以表述为：

非稳态渗流过程的出气端边界条件可以表述为：

U(H,t)＝P_atm+b (6)

非稳态渗流过程的进气端边界条件，可以通过物质守恒原理获得：

其中，P_atm为大气压，可以取为1×10⁵Pa。

步骤(4)中获得控制方程以及上述初始状态和边界条件方程，组成一个二维非线性抛物线型偏微分方程定解问题；对于此类非线性问题，采用有限体积法中的全隐式Picard-Newton迭代格式进行数值求解；Picard-Newton迭代格式具有良好的二阶收敛特性；为构建Picard-Newton迭代格式，对控制方程左边时间项进行向后欧拉时间离散；对控制方程右边的k_∞U/μ项进行Taylor级数展开，继而完成Newton线性化，获得如下线性迭代格式：

式中，j为当前时间步标号，可以取为1,2,…,n-1；n为总时间步数；j＝0对应于初始状态；(s)和(s+1)分别代表第s次和第s+1次迭代计算结果，s可以取为0,1,2,…；Δt为时间步长(s)；

采用全隐Picard空间离散格式，则上式最终可以变换为：

式中：i为当前空间网格点编号，可以取为1,2,…,m-1；i＝0为进气端；i＝m为出气端；为第i个空间网格点在j-1时间步上的最终解，以下表达式中未带迭代计算次数(s)的U值均为最终解。

10.根据权利要求9所述的一种多孔材料的联合测试方法，其特征在于：为启动公式(9)的计算，需要定义所有网格点的初始压力值；

代表初始条件的方程可以离散为：

公式(9)形成含有m+1个未知量的m-1个方程，需要带入上下边界的两个边界条件方程，才能完成计算；代表边界条件的方程可以离散为：

其中

同时，需要假定第j时间步的初始值为第j-1时间步的最终解，即：

根据上述离散的控制方程(公式(9))、初始条件方程(公式(10)和(11))和边界条件方程(公式(12)和(13))，可以形成m+1个线性方程，继而可以求解此线性方程组m+1个未知量以此类推可以求解获得当两次迭代计算结果差值的最大值小于精度要求时，即认为获得第j时间步的最终解；其中，ξ为气压最小分辨率，取为10^-8P₀。

11.根据权利要求10所述的一种多孔材料的联合测试方法，其特征在于：所述步骤(6)的具体步骤为：

对于稳定渗流，步骤(4)控制方程左侧的时间域偏微分项为0，控制方程具有解析解：

将上式以及步骤(2)获得的P₁压力条件下的流量Q带入步骤(4所示的渗流达西定律，可以获得绝对渗透率k_∞与滑移因子b的换算关系：

其中，H为试样高度(m)。

12.根据权利要求11所述的一种多孔材料的联合测试方法，其特征在于：所述步骤(7)的具体步骤为：

假定滑移因子为b、有效孔隙率为φ，并根据公式(16)计算绝对渗透率k_∞；利用步骤(5)中公式(9)、(10)、(11)、(12)和(13)形成线性方程组，计算获得的第j时间步的储气空间压力为采用最小二乘原理设定目标函数如下：

式中，为步骤(2)获得的P(0,t)v.s.t曲线在j时间步处的实测压力值；采用不需要直接计算梯度的修正Powell方向加速法进行计算上述目标函数的最小值，继而获得绝对渗透率k_∞、滑移因子b和有效孔隙率φ最终结果。

13.根据权利要求12所述的一种多孔材料的联合测试方法，其特征在于：为求解上述目标函数最小值及其对应的绝对渗透率k_∞、滑移因子b和有效孔隙率φ，采用不需要直接计算梯度的修正Powell方向加速法进行计算；主要计算步骤为：

步骤7.1：

定义算法起始点v₀ ⁽¹⁾＝[log(b_r),log(φ_r)]、初始搜索方向s₁ ⁽⁰⁾＝[1,0]和s₂ ⁽⁰⁾＝[0,1]、收敛精度ε＝10^-8；其中，b_r和φ_r分别为设定的滑移因子和有效孔隙率的初始值，根据测试材料的经验值进行确定，对于花岗岩可以分别取为1.0×10^-19m²和0.005；

步骤7.2：

对于第j轮计算(j＝1,2,…)，沿方向s₁ ^(j)和s₂ ^(j)逐次进行一维搜索，即

步骤7.3：

分别进行以下计算

判断下述两个条件是否同时满足

f₃＜f₁ (27)

步骤7.4：

当上述两个条件同时满足时，进行以下线性搜索

令

如果停止计算，最终参数值为

否则，重复步骤7.2和7.3进行第j+1轮计算；

步骤7.5：

当步骤7.3中两个条件不能同时满足时，则取

如果停止计算，最终参数值为

否则，重复步骤7.2和7.3进行第j+1轮计算；

步骤7.6：

根据步骤7.4或7.5获得的滑移因子b，利用公式(16)计算绝对渗透率k_∞。