EI、Scopus 收录
中文核心期刊

低渗透煤层气藏中气-水两相不稳定渗流动态分析

刘文超, 刘曰武

刘文超, 刘曰武. 低渗透煤层气藏中气-水两相不稳定渗流动态分析[J]. 力学学报, 2017, 49(4): 828-835. DOI: 10.6052/0459-1879-17-034
引用本文: 刘文超, 刘曰武. 低渗透煤层气藏中气-水两相不稳定渗流动态分析[J]. 力学学报, 2017, 49(4): 828-835. DOI: 10.6052/0459-1879-17-034
Liu Wenchao, Liu Yuewu. DYNAMIC ANALYSIS ON GAS-WATER TWO-PHASE UNSTEADY SEEPAGE FLOW IN LOW-PERMEABLE COALBED GAS RESERVOIRS[J]. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(4): 828-835. DOI: 10.6052/0459-1879-17-034
Citation: Liu Wenchao, Liu Yuewu. DYNAMIC ANALYSIS ON GAS-WATER TWO-PHASE UNSTEADY SEEPAGE FLOW IN LOW-PERMEABLE COALBED GAS RESERVOIRS[J]. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(4): 828-835. DOI: 10.6052/0459-1879-17-034
刘文超, 刘曰武. 低渗透煤层气藏中气-水两相不稳定渗流动态分析[J]. 力学学报, 2017, 49(4): 828-835. CSTR: 32045.14.0459-1879-17-034
引用本文: 刘文超, 刘曰武. 低渗透煤层气藏中气-水两相不稳定渗流动态分析[J]. 力学学报, 2017, 49(4): 828-835. CSTR: 32045.14.0459-1879-17-034
Liu Wenchao, Liu Yuewu. DYNAMIC ANALYSIS ON GAS-WATER TWO-PHASE UNSTEADY SEEPAGE FLOW IN LOW-PERMEABLE COALBED GAS RESERVOIRS[J]. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(4): 828-835. CSTR: 32045.14.0459-1879-17-034
Citation: Liu Wenchao, Liu Yuewu. DYNAMIC ANALYSIS ON GAS-WATER TWO-PHASE UNSTEADY SEEPAGE FLOW IN LOW-PERMEABLE COALBED GAS RESERVOIRS[J]. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(4): 828-835. CSTR: 32045.14.0459-1879-17-034

低渗透煤层气藏中气-水两相不稳定渗流动态分析

基金项目: 

国家自然科学基金 51404232

中国博士后基金 2014M561074

国家科技重大专项 2011ZX05038003

详细信息
    通讯作者:

    2) 刘文超, 讲师, 主要研究方向:渗流力学.E-mail:wcliu 2008@126.com

    3) 刘曰武, 研究员, 主要研究方向:渗流力学.E-mail:lywu@imech.ac.cn

  • 中图分类号: O357.3

DYNAMIC ANALYSIS ON GAS-WATER TWO-PHASE UNSTEADY SEEPAGE FLOW IN LOW-PERMEABLE COALBED GAS RESERVOIRS

  • 摘要: 针对低渗透煤层渗流问题,考虑启动压力梯度及其引起的动边界和动边界内吸附气解吸作用的渗流模型研究目前仅限于单相流,而更符合实际的气-水两相渗流动边界模型未见报道.本文综合考虑了煤层吸附气的解吸作用、气-水两相渗流、非达西渗流、地层应力敏感等影响因素,进行了低渗透煤层的气-水两相渗流模型研究.采用了试井技术中的“分相处理”方法,修正了两相渗流的综合压缩系数和流度,并基于含气饱和度呈线性递减分布的假设,建立了煤层气藏的气-水两相渗流耦合模型.该数学模型不仅可以描述由于低渗透煤层中渗流存在启动压力梯度而产生的可表征煤层有效动用范围随时间变化的移动边界,还可以描述煤层有效动用范围内吸附气的解吸现象以及吸附气解吸作用所引起的煤层含气饱和度的上升;为了提高模型精度,控制方程还保留了二次压力梯度项.采用了稳定的全隐式有限差分方法进行了模型的数值求解,并验证了数值计算方法的正确性,获得了模型关于瞬时井底压力与压力导数响应的双对数特征曲线,由此分析了各渗流参数的敏感性影响.本文研究结果可为低渗透煤层气藏开发的气-水两相流试井技术提供渗流力学的理论基础.
    Abstract: At present, the research on the seepage flow models in low-permeable coalbeds is only limited to the single phase flow cases, which simultaneously considerate the existence of threshold pressure gradient in the seepage flow process, its produced moving boundary and the desorption function of adsorbed gas inside the moving boundary; however, the research on the gas and water two-phase seepage models with moving boundaries, which are more consistent with the actual situations, has not been reported yet. In comprehensive consideration of these influential factors including the desorption function of the adsorbed gas in coalbeds, gas-water two-phase seepage flow, the non-Darcy seepage flow characteristics in the low-permeable formations, the stress-sensitive effect of the formation, etc., modeling the gas-water two-phase seepage flow in low-permeable coalbeds is studied in this paper. According to the "phase separation" method involved in the well testing technology, both the comprehensive compressibility coefficient and the fluidity are modified for the two-phase seepage flow problem; and then based on the assumption of the linear spatial distribution of the gas saturation, a coupled model of gas-water two-phase seepage flow in low-permeable coalbeds is built. The mathematical model can not only depict the moving boundary that represents the change of the effectively disturbed coalbed area with time due to the existence of threshold pressure gradient in the seepage flow process in low-permeable coalbeds, but also can depict the desorption phenomena of the adsorbed gas in the effectively disturbed coalbed area, and the increase of the gas saturation in coalbeds caused by the desorption function of the adsorbed gas; furthermore, in order to improve the accuracy of the model, the quadratic pressure gradient term is retained in the governing equation of the model. A fully implicit finite difference method is adopted to numerically solve the model, and the correctness of the numerical method is also verified. Eventually, the log-log type curves regarding the transient wellbore pressure response and its derivative are obtained from the model, and then the sensitive-effects of some seepage flow parameters can be analyzed. The presented research results in the paper can provide theoretical foundations of seepage flow mechanics for the well testing technology for the gas-water two-phase seepage flow in the development of low-permeable coalbed gas reservoirs.
  • 与常规油气藏相比,煤层气藏有着截然不同的能源储集方式、开采方式及渗流机制:① 煤层气藏中存在大量吸附气[1-4],储集于煤层基质中;② 煤层气开采通常是一个排水降压过程[1-2],生产过程中当地层压力下降至临界解吸压力后,吸附气会解吸并作为物质供给,因而也会涉及到气-水两相流;③ 煤层气藏开采过程中,由于低渗透煤层孔隙狭小、表面积大,流-固界面作用强烈,以及低渗透煤层的非均质性、孔隙边界层以及气-水两相流毛管阻力的影响[4-12],低渗煤层中的两相渗流还存在启动压力梯度[5-6],会导致煤层气藏有效动用范围随时间而变化[13-19], 即:有效动用范围仅存在于压力梯度大于启动压力梯度的动边界以内,而动边界以外不发生流动, 渗流模型计算研究表明[13-19], 考虑启动压力梯度的非达西渗流模型应考虑动边界的影响; ④ 随着地层压力的不断降低,应力敏感的岩石孔隙会发生变形,导致煤层渗透能力发生变化[20].

    目前,对于同时考虑低渗透煤层渗流存在启动压力梯度及其引起的动边界和动边界内吸附气解吸作用的煤层气渗流模型研究仅限于单相渗流动边界模型[21-23],而考虑吸附气解吸作用的气-水两相渗流动边界模型研究仍未见相关报道.本文将综合考虑低渗透煤层吸附气的解吸作用、气-水两相渗流、低渗煤层的非达西渗流特征、地层应力敏感等影响因素,在以往考虑动边界的低渗透煤层单相(包括单相水和单相气)渗流模型研究[21-26]基础上,结合试井解释中的``分相处理''方法[26-27],并通过引入稳定解吸与不稳定解吸的两个源项[22-26]来表达煤层吸附气的解吸作用,进一步建立了一个符合低渗透煤层气排采特征的气-水两相渗流耦合模型.本文研究可为低渗透煤层气藏开发所涉及的气-水两相流试井技术提供渗流力学的理论基础.

    煤层气生产时,生产井周围地层压力最先开始低于吸附气的临界解吸压力,煤层气发生解吸后,不断流动至裂隙系统,气-水混合后发生由裂隙系统向井筒的气-水两相渗流,煤层气藏的气-水两相渗流特征如图 1所示.由于低渗煤层存在启动压力梯度,煤层渗流存在动边界[13, 21-23],井筒至动边界之间的渗流区域属于煤层压降区域,地层压力梯度大于启动压力梯度流体发生流动,表征煤层的有效动用范围;而动边界外为煤层未动用区域,不存在压力降.假设临界解吸压力与原始地层压力大致相等,则可近似认为煤层压降区域内都发生解吸;随着时间的延长,煤层吸附气会不断解吸到煤层中,动边界以内煤层含气饱和度不断增大[26].

    图  1  低渗透煤层气藏气-水两相渗流模型示意图
    Figure  1.  The schematics of the gas-water two-phase seepage flow model in low-permeable coalbed gas reservoirs

    模型假设条件为[21-23]

    (1) 煤层气储层为均质、各向同性的无限大地层;忽略重力及毛管力影响,不考虑其他的物理化学因素;不考虑温度变化对流动的影响.

    (2) 低渗煤层中渗流规律遵循考虑启动压力梯度的非达西渗流运动方程[5-6, 13]

    $$ \upsilon = \left\{ \begin{array}{l} - \frac{k}{\mu } \cdot \frac{{\partial p}}{{\partial r}} \cdot \left( {1 - \lambda /\left| {\frac{{\partial p}}{{\partial r}}} \right|} \right){\mkern 1mu} \left| {\frac{{\partial p}}{{\partial r}}} \right| > \lambda \\ 0, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0\left| {\frac{{\partial p}}{{\partial r}}} \right|\lambda \end{array} \right. $$ (1)

    式中, $\lambda $ 为气-水两相渗流的启动压力梯度; $r$ 为径向距离; $p$ 为地层压力; $k$ 为岩石渗透率; $\upsilon $ 为渗流速度; $\mu $ 为流体黏度.

    (3) 假定压敏效应引起的煤层渗透率变化为关于地层压力的指数函数[20-21]

    $$ k = k_{0} {\rm e}^{ - \gamma \left( {p_0 - p} \right)} $$ (2)

    (4) 基质表面解吸的煤层气直接进入裂隙,通过引入临界解吸压力和解吸系数,在常规的控制方程中添加恒定的源项(稳定解吸)和非恒定的源项(不稳定解吸)来表达煤层气的解吸作用[22-26],建立一种考虑煤层气解吸的均匀介质模型.

    (5) 假定煤层气临界解吸压力 $p_{\rm c}$ 与初始地层压力大致相等,煤层压降区域内都发生解吸.

    (6) 煤层气解吸后流向井筒,越靠近煤层气井,含气饱和度越高,而动边界以外地层压力保持为初始压力,为未解吸区域;参照文献[26]中提出的两区线性分布模型,假定从井筒到动边界 $\delta $ 处(对应启动压力梯度的动边界位置)的范围内煤层气藏的含气饱和度服从线性递减分布[26](图 2);假定井底含气饱和度 $S_{g\_{\rm well}}$ 在一定的生产时间内保持不变;动边界处含气饱和度为初始含气饱和度 $S_{g\_{\rm ini}}$ .

    图  2  含气饱和度分布随时间变化示意图
    Figure  2.  The schematics of the change of the gas saturation distribution with time

    根据试井解释技术中分相处理方法[26-27],建立了煤层气藏的气-水两相渗流模型.在单相渗流相关公式基础[21-23]上,基于气-水两相渗流方程,从综合弹性压缩系数和流度两个关键参数出发,进行了气-水两相渗流公式的修正.

    修正后的气-水两相渗流模型的综合压缩系数可表示为[26-27]

    $$ C_{\rm t} = C_{\rm f} + S_{\rm g} \cdot C_{\rm g} + \left( {1 - S_{\rm g} } \right) \cdot C_{\rm w} $$ (3)

    式中, $S_{\rm g}$ 为煤层的含气饱和度; $C_{\rm g}$ 为气体的等温压缩系数; $C_{\rm w}$ 为水的压缩系数.

    修正后的气-水两相渗流模型的流度 $\lambda_{\rm t}$ 可表示为[26-27]

    $$ \lambda _{\rm t} = \dfrac{k \cdot k_{{\rm rg}} \left( {S_{\rm g} } \right)}{\mu _{\rm g} } + \dfrac{k \cdot k_{{\rm rw}} \left( {S_{\rm g} } \right)}{\mu _{\rm w} } $$ (4)

    式中, $\mu_{\rm g}$ 为气相黏度; $\mu_{\rm w}$ 为水相黏度; $k_{\rm rg}$ 为气相相对渗透率; $k_{\rm rw}$ 为水相相对渗透率.

    低渗透煤层有效动用区域内存在吸附气的稳定解吸与不稳定解吸,并通过在控制方程中分别添加恒定的源项 $\alpha_{\rm 1D}$ 和非恒定的源项 $\alpha_{\rm 2D}\cdot (p_{\rm D}-p_{\rm cD})$ 来表示;其中,稳定解吸代表煤层吸附气的稳定解吸作用,而不稳定解吸与压力的降低程度有关[2, 22-24];由此推导了低渗透煤层气藏的气-水两相渗流耦合模型.

    模型的控制方程为

    $$ \begin{array}{l} \frac{1}{{{r_{\rm{D}}}}} \cdot \frac{\partial }{{\partial {r_{\rm{D}}}}}\left[{{{\rm{e}}^{-{\gamma _{\rm{D}}} \cdot {p_{\rm{D}}}}} \cdot {r_{\rm{D}}} \cdot \left( {\frac{{\partial {p_{\rm{D}}}}}{{\partial {r_{\rm{D}}}}} + {\lambda _{\rm{D}}}} \right)} \right] + \\ \left[{{\alpha _{{\rm{1D}}}} + {\alpha _{{\rm{2D}}}}\left( {{p_{\rm{D}}}-{p_{{\rm{cD}}}}} \right)} \right] = \\ \frac{{{C_{{\rm{tD}}}}\left( {{S_{\rm{g}}}} \right)}}{{{\lambda _{{\rm{tD}}}}\left( {{S_{\rm{g}}}} \right)}} \cdot \frac{{\partial {p_{\rm{D}}}\left( {{r_{\rm{D}}}, {t_{\rm{D}}}} \right)}}{{\partial {t_{\rm{D}}}}}{\mkern 1mu}, \;1 \le {r_{\rm{D}}} \le \delta \left( {{t_{\rm{D}}}} \right) \end{array} $$ (5)

    其中

    $$ \begin{array}{*{35}{l}} {{\lambda }_{\text{tD}}}\left( {{S}_{\text{g}}} \right)=\frac{{{k}_{\text{rg}}}\left( {{S}_{\text{g}}} \right)/{{\mu }_{\text{g}}}+{{k}_{\text{rw}}}\left( {{S}_{\text{g}}} \right)/{{\mu }_{\text{w}}}}{{{k}_{\text{rw}}}\left( {{S}_{\text{g}\_\text{initial}}} \right)/{{\mu }_{\text{w}}}+{{k}_{\text{rg}}}\left( {{S}_{\text{g}\_\text{initial}}} \right)/{{\mu }_{\text{g}}}} \\ {{C}_{\text{tD}}}\left( {{S}_{\text{g}}} \right)=\frac{{{C}_{\text{f}}}+{{S}_{\text{g}}}\cdot {{C}_{\text{g}}}+\left( 1-{{S}_{\text{g}}} \right)\cdot {{C}_{\text{w}}}}{{{C}_{\text{f}}}+{{S}_{\text{g}\_\text{initial}}}\cdot {{C}_{\text{g}}}+\left( 1-{{S}_{\text{g}\_\text{initial}}} \right)\cdot {{C}_{\text{w}}}} \\ {{S}_{\text{g}}}\left( {{r}_{\text{D}}},{{t}_{\text{D}}} \right)=\left\{ \begin{array}{*{35}{l}} \begin{array}{*{35}{l}} {{S}_{\text{g}\_\text{well}}}-\frac{{{S}_{\text{g}\_\text{well}}}-{{S}_{\text{g}\_\text{ini}}}}{\delta \left( {{t}_{\text{D}}} \right)-1}\cdot \left( {{r}_{\text{D}}}-1 \right), \\ \ 1\le {{r}_{\text{D}}}\le \delta ,\ \\ \ {{S}_{\text{g}\_\text{ini}}},{{r}_{\text{D}}}>\delta \\ \end{array} \\ \end{array} \right. \\ \end{array} $$

    为了提高模型精度,式(5) 保留了二次压力梯度项 $\gamma_{\rm D}\cdot (\partial p_{\rm D}/\partial r_{\rm D})^{2}$ [15, 21, 28].建模过程中将依赖于地层压力的指数形式, 即式(2), 代入质量守恒方程时会产生该项.

    初始条件为

    $$ \left. {p_{\rm D} } \right|_{t_{\rm D} = 0} = 0 $$ (6)
    $$ \left. {S_{\rm g} } \right|_{t_{\rm D} = 0} = S_{{\rm g\_initial}} $$ (7)

    考虑井筒储集和表皮效应的定流量内边界条件为[23, 26]

    $$ C_{\rm D} \cdot \dfrac{{\rm d}p_{{\rm WD}} }{{\rm d}t_{\rm D} } - \left. {{\rm e}^{ - \gamma _{\rm D} \cdot p_{\rm D} } \cdot A_{\rm \lambda } \cdot \left( {\dfrac{\partial p_{\rm D} }{\partial r_{\rm D} } + \lambda _{\rm D} } \right)} \right|_{r_{\rm D} = 1} = 1 $$ (8)
    $$ P_{{\rm WD}} = \left. {\left[{p_{\rm D}-S \cdot {\rm e}^{-\gamma _{\rm D} \cdot p_{\rm D} } \cdot A_{\rm \lambda } \cdot \left( {\dfrac{\partial p_{\rm D} }{\partial r_{\rm D} } + \lambda _{\rm D} } \right)} \right]} \right|_{r_{\rm D} = 1} $$ (9)

    其中

    $$ A_{\rm \lambda } = \dfrac{k_{{\rm rw}} \left( {S_{{\rm g\_well}} } \right) / \mu _{\rm w} + k_{{\rm rg}} \left( {S_{{\rm g\_well}} } \right) / \mu _{\rm g} }{k_{{\rm rw}} \left( {S_{{\rm g\_initial}} } \right) / \mu _{\rm w} + k_{{\rm rg}} \left( {S_{{\rm g\_initial}} } \right) / \mu _{\rm g} } $$

    表征煤层有效动用范围随时间而变化的动边界条件为[21-23]

    $$ \left. {\dfrac{\partial p_{\rm D} }{\partial r_{\rm D} }} \right|_{r_{\rm D} = \delta \left( {t_{\rm D} } \right)} = - \lambda _{\rm D} $$ (10)
    $$ \left. {p_{\rm D} } \right|_{r_{\rm D} = \delta \left( {t_{\rm D} } \right)} = 0 $$ (11)
    $$ \left. {S_{\rm g} } \right|_{r_{\rm D} = \delta \left( {t_{\rm D} } \right)} = S_{{\rm g\_initial}} $$ (12)

    式中, $p_{\rm D}$ 为无因次地层压力; $p_{\rm cD}$ 为无因次临界解吸压力; $r_{\rm D}$ 为无因次径向距离; $t_{\rm D}$ 为无因次时间; $\gamma_{\rm D}$ 为表征岩石压敏效应对煤层渗透率影响的无因次渗透率模量; $\lambda_{\rm D}$ 为无因次启动压力梯度; $\delta $ 为无因次动边界距离; $\alpha_{\rm 1D}$ 为无因次稳定解吸系数[22-26]; $\alpha_{\rm 2D}$ 为无因次不稳定解吸系数[22-26]; $p_{\rm WD}$ 为无因次井底压力; $C_{\rm D}$ 无因次井筒储集系数; $S$ 为表皮因子.

    无量纲变量定义为

    $$ \begin{array}{l} \tilde \mu = {\left[{{k_{{\rm{rw}}}}\left( {{S_{{\rm{g}}\_{\rm{initial}}}}} \right)/{\mu _{\rm{w}}} + {k_{{\rm{rg}}}}\left( {{S_{{\rm{g}}\_{\rm{initial}}}}} \right)/{\mu _{\rm{g}}}} \right]^{ - 1}}\\ {p_{{\rm{WD}}}} = \frac{{2{\rm{\pi }}{k_0}h}}{{q\tilde \mu }}\left( {{p_0} - {p_{\rm{w}}}} \right){\mkern 1mu}, \;\;{\alpha _{{\rm{1D}}}} = - \frac{{2{\rm{\pi }} \cdot h \cdot r_{\rm{w}}^2}}{q} \cdot {\alpha _1}\\ {\alpha _{{\rm{2D}}}} = - \frac{{\mu \cdot r_{\rm{w}}^2}}{{{k_{\rm{0}}}}} \cdot {\alpha _2}{\mkern 1mu}, \;\;{\lambda _{\rm{D}}} = \frac{{2{\rm{\pi }}{k_{\rm{0}}}h{r_{\rm{w}}}\lambda }}{{q\tilde \mu }}{\mkern 1mu}, \;\;{\gamma _{\rm{D}}} = \frac{{q\tilde \mu \gamma }}{{2{\rm{\pi }}{k_0}h}}\\ {r_{\rm{D}}} = r/{r_{\rm{w}}}{\mkern 1mu}, \;{t_{\rm{D}}} = \frac{{{k_0}}}{{{\phi _0}{C_{\rm{t}}}\tilde \mu r_{\rm{w}}^2}}t{\mkern 1mu}, \;\;{C_{\rm{D}}} = \frac{C}{{2{\rm{\pi }}{\phi _0}{C_{\rm{t}}}hr_{\rm{w}}^2}}\\ {p_{\rm{D}}}\left( {{r_{\rm{D}}}, {t_{\rm{D}}}} \right) = \frac{{2{\rm{\pi }}{k_0}h}}{{q\tilde \mu }}\left[{{p_0}-p\left( {r, t} \right)} \right]\\ {p_{{\rm{cD}}}} = \frac{{2{\rm{\pi }}{k_0}h}}{{q\tilde \mu }}\left( {{p_0} - {p_{\rm{c}}}} \right) \end{array} $$

    其中, $\alpha_{1}$ 和 $\alpha_{2}$ 分别为煤层气稳定解吸系数和不稳定解吸系数[22-26]; $\tilde \mu$ 为两相流体的等效黏度; $h$ 为煤层厚度; $r_{w}$ 为井筒半径; $\phi_{0}$ 为煤层的初始孔隙度; $C$ 为井筒储集系数; $q$ 为地层条件下的流体流量; $p_{\rm w}$ 为井底压力; $p_{\rm c}$ 为临界解吸压力[22-26].上述所有有量纲的参数单位均为统一单位制.

    由于数学模型存在动边界,渗流区域并非固定,而是随着时间逐渐向外扩展;在数值求解时,为克服瞬时流动区域网格划分的难题,引入如下空间坐标变换[21, 28-29]

    $$ y_{\rm D} \left( {r_{\rm D}, t_{\rm D} } \right) = \dfrac{r_{\rm D} - 1}{\delta \left( {t_{\rm D} } \right) - 1}\, , \ \ 1 \leqslant r_{\rm D} \leqslant \delta \left( {t_{\rm D} } \right) $$ (13)

    式(13) 经坐标变换后,具有动边界的数学模型可转化为具有定边界的数学模型,具体的公式推导方法可参照文献[21].本文采用全隐式的有限差分方法对转化后的非线性数学模型进行数值求解.建立差分方程时,一阶导数采用一阶向前差分,二阶导数采用二阶中心差分.最后,采用Newton-Raphson迭代方法[21]对差分方程进行数值求解.

    当 $C_{\rm tD}/\lambda_{\rm tD}$ , $A_{\lambda }=1.0$ 时,模型可退化为单相情形, 数值结果可与现有的利用其他方法所求得的模型计算结果[30]进行对比,以验证本文数值计算方法的正确性.图 3为 $\lambda _{\rm D}=0.1$ 时,在双对坐标中关于无因次瞬时井底压力及压力导数的模型计算结果对比.由图 3可以看出,本文与文献[30]的数值结果吻合较好,由此验证了本文计算方法的正确性.

    图  3  数值计算方法的正确性验证
    Figure  3.  The verification of the correctness of the numerical method

    根据煤层气生产实际情况,煤层气-水两相渗流模型的无因次的基础参数数据为: $\lambda_{\rm D}=0.2$ , $\alpha_{\rm 1D}=-0.000\, 01$ , $\alpha_{\rm 2D}=-0.000\, 01$ , $P_{\rm cD} =0$ , $\gamma_{\rm D}=0.001$ , $C_{\rm D} =20.0$ , $S =10.0$ ,束缚水饱和度 $S_{\rm wc} =0.15$ , $S_{\rm g\_ini} =0.1$ , $S_{\rm g\_well} =0.65$ , $C_{\rm g}/C_{\rm w} =25$ , $C_{\rm f}/C_{\rm w} =0.5$ , $\mu _{\rm w}/\mu_{\rm g}=100$ ;水相和气相的相对渗透率计算公式[26]

    $$ \begin{array}{l} {k_{{\rm{rg}}}} = \left( {1 - \frac{{1 - {S_{\rm{g}}} - {S_{{\rm{wc}}}}}}{{1 - {S_{{\rm{wc}}}}}}} \right) \cdot \\ \sqrt {1 - {{\left( {\frac{{1 - {S_{\rm{g}}} - {S_{{\rm{wc}}}}}}{{1 - {S_{{\rm{wc}}}}}}} \right)}^{{\textstyle{1 \over 4}}}} \cdot {{\left( {1 - {S_{\rm{g}}} - {S_{{\rm{wc}}}}} \right)}^{{\textstyle{1 \over 2}}}}} \\ {k_{{\rm{rw}}}} = {\left( {1 - {S_{\rm{g}}}} \right)^3} \cdot {\left( {\frac{{1 - {S_{\rm{g}}} - {S_{{\rm{wc}}}}}}{{1 - {S_{{\rm{wc}}}}}}} \right)^{{\textstyle{3 \over 2}}}} \end{array} $$

    采用上述数据作为模型输入参数进行了煤层气-水两相渗流耦合模型的计算,并对模型的数值计算结果进行了分析.

    由于考虑了煤层压敏效应所引起的煤层渗透率的指数变化,在推导煤层气渗流模型时,控制方程式(5) 会产生二次压力梯度项[15, 21, 28] ${{\gamma }_{\text{D}}}{{(\partial {{P}_{\text{D}}}/\partial {{r}_{\text{D}}})}^{2}}$ , 其系数为无因次渗透率模量 $\gamma_{\rm D}$ .图 4为忽略二次压力梯度项对无因次瞬时井底压力及压力导数的影响曲线.由图 4可以看出,当 $\gamma_{\rm D}=0.015$ 时,保留了二次压力梯度项才可反映出无因次渗透率模量 $\gamma_{\rm D}$ 对无因次瞬时井底压力及压力导数的影响;无因次时间越长,影响越明显;然而,如果忽略二次压力梯度项, $\gamma _{\rm D}=0.015$ 与 $\gamma_{\rm D}=0$ 分别对应的无因次瞬时井底压力与压力导数曲线将重合在一起.

    图  4  忽略二次压力梯度项的影响
    Figure  4.  The effect of neglecting the quadratic pressure gradient term

    图 5为 $\gamma_{\rm D} =0.015$ 时,忽略二次压力梯度项所引起的瞬时井底压力及压力导数的绝对误差与相对误差.由图 5可以看出,忽略二次压力梯度项对瞬时井底压力与压力导数计算的影响显著,造成压力导数计算的相对误差更大;无因次时间越长,瞬时井底压力与压力导数的相对误差也会越大,压力导数计算的最大相对误差超过了20%.因此,该低渗透煤层气-水两相渗流耦合模型应考虑二次压力梯度项的影响.

    图  5  忽略二次压力梯度项所引起的绝对误差与相对误差
    Figure  5.  The absolute error and the relative error caused by neglecting the quadratic pressure gradient term

    图 6为不同时刻下,关于无因次综合压缩系数 $C_{\rm tD}$ 与无因次流度 $\lambda_{\rm tD}$ 比值的空间分布曲线.由图 6可以看出,在煤层的有效动用区域内, $C_{\rm tD}$ 与 $\lambda_{\rm tD}$ 的比值小于1.0;由曲线的变化特征还可看出,径向距离 $r_{\rm D}$ 越小, $C_{\rm tD}$ 与 $\lambda _{\rm tD}$ 的比值也越小,气-水两相渗流特征越明显;随时间的延长, $C_{\rm tD}$ 与 $\lambda_{\rm tD}$ 比值也越小,气-水两相渗流特征也越明显,这是由于随着地层压力持续降低,煤层吸附气不断解吸至含水的煤层中;然而, $C_{\rm tD}$ 与 $\lambda _{\rm tD}$ 比值变化幅度会越来越小.

    图  6  无因次综合压缩系数 $C_{\rm tD}$ 与无因次流度 $\lambda_{\rm tD}$ 比值的变化
    Figure  6.  The change of the ratio of the dimensionless comprehensive compressibility $C_{\rm tD}$ to the dimensionless fluidity $\lambda_{\rm tD}$

    图 7为无因次稳定解吸系数 $\alpha_{\rm 1D}$ 与无因次不稳定解吸系数 $\alpha_{\rm 2D}$ 对无因次瞬时井底压力及压力导数双对数特征曲线影响.由图 7可以看出,煤层气解吸作用会使无因次瞬时井底压力导数曲线的后期下掉,无因次瞬时井底压力降低.这是由于在煤层气稳定解吸和不稳定解吸的作用下,煤层气不断解吸出来,及时补充到煤层气藏中,减缓了地层压力的降落,导致无因次瞬时井底压力和无因次井底压力导数的下降.煤层气解吸作用的这一影响特征已由实际煤层气井试井数据制作出的井底压力及压力导数的双对数特征曲线所证实[31].

    图  7  无因次吸附气解吸系数的影响
    Figure  7.  The effect of the dimensionless coefficient of the desorption of adsorbed gas

    图 8为井底含气饱和度 $S_{\rm g\_well}$ 对无因次瞬时井底压力及压力导数双对数特征曲线的影响.由图 8可以看出,井底含气饱和度 $S_{\rm g\_well}$ 越大,无因次瞬时井底压力及井底压力导数越小;无因次时间越长,井底含气饱和度 $S_{\rm g\_well}$ 的影响会越明显.由煤层含气饱和度的表达式可知,井底含气饱和度 $S_{\rm g\_well}$ 越大,煤层压降区域内含气饱和度也越大,表征流体流动能力无因次流度值也会增大;综合压缩系数也会增大,但不足以抵消煤层中流体流动能力提升对两相渗流的影响,因而延缓了地层压力的降落,地层压力降落的速度也变慢.

    图  8  含气饱和度的影响
    Figure  8.  The effect of the gas saturation

    图 9为无因次启动压力梯度 $\lambda_{\rm D}$ 对无因次瞬时井底压力及压力导数双对数特征曲线的影响.图 9表明, 无因次启动压力梯度越大,无因次井底压力和压力导数也越大,双对数曲线向上抬升,而无因次启动压力梯度影响的敏感性会减弱.这是由于启动压力梯度的存在会造成消耗更多的地层能量,加剧了地层压力的降落.

    图  9  无因次启动压力梯度的影响
    Figure  9.  The effect of the dimensionless threshold pressure gradient

    (1) 采用试井解释中的分相处理方法,建立了低渗透煤层气藏的气-水两相渗流耦合模型.该模型不仅可以描述煤层气生产时由于启动压力梯度存在而产生的可表征煤层有效动用范围随时间变化的移动边界,而且还可以考虑有效动用区域内煤层气的稳定解吸与不稳定解吸作用,以及吸附气解吸所引起的煤层含气饱和度的变化.

    (2) 采用了稳定的全隐式有限差分方法进行了模型的数值求解,并进行了数值计算方法的正确性验证.

    (3) 模型计算结果表明:① 保留二次压力梯度项才可反映出无因次渗透率模量对无因次瞬时井底压力与压力导数的影响;忽略二次压力梯度项会产生较大的相对误差. ② 径向距离越小, $C_{\rm tD}$ 与 $\lambda_{\rm tD}$ 比值也越小,气-水两相渗流的特征更明显. ③ 煤层气解吸作用会使无因次瞬时井底压力导数曲线的后期下掉. ④ 煤层压降区域内含气饱和度越大,流体在煤层中流动能力会越强,因而会延缓地层压力的降落. ⑤ 启动压力梯度越大,无因次瞬时井底压力与压力导数的双对数曲线向上抬升,然而其影响的敏感性会减弱.

  • 图  1   低渗透煤层气藏气-水两相渗流模型示意图

    Figure  1.   The schematics of the gas-water two-phase seepage flow model in low-permeable coalbed gas reservoirs

    图  2   含气饱和度分布随时间变化示意图

    Figure  2.   The schematics of the change of the gas saturation distribution with time

    图  3   数值计算方法的正确性验证

    Figure  3.   The verification of the correctness of the numerical method

    图  4   忽略二次压力梯度项的影响

    Figure  4.   The effect of neglecting the quadratic pressure gradient term

    图  5   忽略二次压力梯度项所引起的绝对误差与相对误差

    Figure  5.   The absolute error and the relative error caused by neglecting the quadratic pressure gradient term

    图  6   无因次综合压缩系数 $C_{\rm tD}$ 与无因次流度 $\lambda_{\rm tD}$ 比值的变化

    Figure  6.   The change of the ratio of the dimensionless comprehensive compressibility $C_{\rm tD}$ to the dimensionless fluidity $\lambda_{\rm tD}$

    图  7   无因次吸附气解吸系数的影响

    Figure  7.   The effect of the dimensionless coefficient of the desorption of adsorbed gas

    图  8   含气饱和度的影响

    Figure  8.   The effect of the gas saturation

    图  9   无因次启动压力梯度的影响

    Figure  9.   The effect of the dimensionless threshold pressure gradient

  • [1] 刘曰武, 赵培华, 鹿倩等.煤层气与常规天然气测试技术的异同.油气井测试, 2010, 19(6): 6-11 http://www.cnki.com.cn/Article/CJFDTOTAL-YQJC201006002.htm

    Liu Yuewu, Zhao Peihua, Lu Qian, et al. Differences of well testing between CBM and conventional gas well. Well Testing, 2010, 19(6): 6-11 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-YQJC201006002.htm

    [2] 刘曰武, 苏中良, 方虹斌等.煤层气的解吸/吸附机理研究综述.油气井测试, 2010, 19(6): 37-44 http://www.cnki.com.cn/Article/CJFDTOTAL-YQJC201006007.htm

    Liu Yuewu, Su Zhongliang, Fang Hongbin, et al. Review on CBM desorption/adsorption mechanism. Well Testing, 2010, 19(6): 37-44 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-YQJC201006007.htm

    [3] 郭为, 胡志明, 左罗等.页岩基质解吸-扩散-渗流耦合实验及数学模型.力学学报, 2015, 47(6): 916-922 doi: 10.6052/0459-1879-15-068

    Guo Wei, Hu Zhiming, Zuo Luo, et al. Gas desorption-diffusion-seepage coupled experiment of shale matrix and mathematic model. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(6): 916-922 (in Chinese) doi: 10.6052/0459-1879-15-068

    [4] 李靖, 李相方, 王香增等.页岩黏土孔隙含水饱和度分布及其对甲烷吸附的影响.力学学报, 2016, 48(5): 1217-1228 http://lxxb.cstam.org.cn/CN/article/downloadArticleFile.do?attachType=PDF&id=145828

    Li Jing, Li Xiangfang, Wang Xiangzeng, et al. Effect of water distribution on methane adsorption capacity in shale clay. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(5): 1217-1228 (in Chinese) http://lxxb.cstam.org.cn/CN/article/downloadArticleFile.do?attachType=PDF&id=145828

    [5] 郭红玉, 苏现波.煤储层启动压力梯度的实验测定及意义.天然气工业, 2010, 30(6): 52-54 doi: 10.3787/j.issn.10000976.2010.06.014

    Guo Hongyu, Su Xianbo. An experimental measurement of the threshold pressure gradient of coal reservoirs and its significance. Natural Gas Industry, 2010, 30(6): 52-54 (in Chinese) doi: 10.3787/j.issn.10000976.2010.06.014

    [6] 田巍, 朱维耀, 朱华银等.致密砂岩凝析气藏启动压力梯度.中南大学学报(自然科学版), 2015, 46(9): 3415-3421 doi: 10.11817/j.issn.1672-7207.2015.09.034

    Tian Wei, Zhu Weiyao, Zhu Huayin, et al. Condense gas reservoir threshold pressure gradient for tight sandstone. Journal of Central South University (Science and Technology), 2015, 46(9): 3415-3421 (in Chinese) doi: 10.11817/j.issn.1672-7207.2015.09.034

    [7] 朱维耀, 宋洪庆, 何东博等.含水低渗气藏低速非达西渗流数学模型及产能方程研究.天然气地球科学, 2008, 19(5): 685-689 http://www.cnki.com.cn/Article/CJFDTOTAL-TDKX200805021.htm

    Zhu Weiyao, Song Hongqing, He Dongbo, et al. Low-velocity non-Darcy gas seepage model and productivity equations of low-permeability water-bearing gas reservoirs. Natural Gas Geoscience, 2008, 19(5): 685-689 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-TDKX200805021.htm

    [8] 宋付权, 刘慈群, 吴柏志.启动压力梯度的不稳定快速测量.石油学报, 2001, 22(3): 67-70 doi: 10.7623/syxb200103014

    Song Fuquan, Liu Ciqun, Wu Bozhi. The fast transient measurement of threshold pressure gradient. Acta Petrolei Sinica, 2001, 22(3): 67-70 (in Chinese) doi: 10.7623/syxb200103014

    [9] 李爱芬, 张少辉, 刘敏等.一种测定低渗油藏启动压力的新方法.中国石油大学学报(自然科学版), 2008, 32(1): 68-71 http://www.cnki.com.cn/Article/CJFDTOTAL-SYDX200801019.htm

    Li Aifen, Zhang Shaohui, Liu Min, et al. A new method of measuring starting pressure for low permeability reservoir. Journal of China University of Petroleum, 2008, 32(1): 68-71 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-SYDX200801019.htm

    [10] 李爱芬, 刘敏, 张化强等.低渗透油藏油水两相启动压力梯度变化规律研究.西安石油大学学报(自然科学版), 2010, 25(6): 47-54 http://cdmd.cnki.com.cn/Article/CDMD-10220-2004050749.htm

    Li Aifen, Liu Min, Zhang Huaqiang, et al. Experimental study on threshold pressure gradient of oil-water two-phase seepage in low permeability reservoirs. Journal of Xi'an Shiyou University (Natural Science Edition), 2010, 25(6): 47-54 (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10220-2004050749.htm

    [11] 刘曰武, 丁振华, 何凤珍.确定低渗透油藏启动压力梯度的三种方法.油气井测试, 2002, 11(4): 1-4 http://www.cnki.com.cn/Article/CJFDTOTAL-YQJC200204000.htm

    Liu Yuewu, Ding Zhenhua, He Fengzhen. Three kinds of methods for determining the start-up pressure gradients in low permeability reservoir. Well Testing, 2002, 11(4): 1-4 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-YQJC200204000.htm

    [12] 舒卫兵, 许鹤华, 刘唐伟等.低渗透裂缝性油气藏非稳态窜流因子研究.力学学报, 2014, 46(1): 70-77 doi: 10.6052/0459-1879-13-154

    Shu Weibing, Xu Hehua, Liu Tangwei, et al. Research on nonsteady shape factor of cross flow in low permeability fractured reservoirs. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(1): 70-77 (in Chinese) doi: 10.6052/0459-1879-13-154

    [13]

    Yao J, Liu WC, Chen ZX. Numerical solution of a moving boundary problem of one-dimensional flow in semi-infinite long porous media with threshold pressure gradient. Mathematical Problems in Engineering], 2013, 2013: 384246.

    [14] 宋付权, 刘慈群, 李凡华.低渗透介质含启动压力梯度一维瞬时压力分析.应用数学和力学, 1999, 20(1): 25-32 http://www.cnki.com.cn/Article/CJFDTOTAL-YYSX901.003.htm

    Song Fuquan, Liu Ciqun, Li Fanhua. Transient pressure of percolation through one dimension porous media with threshold pressure gradient. Applied Mathematics and Mechanics, 1999, 20(1): 25-32 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-YYSX901.003.htm

    [15]

    Liu WC, Yao J. Numerical investigations of the effect of nonlinear quadratic pressure gradient term on a moving boundary problem of radial flow in low-permeable reservoirs with threshold pressure gradient. Mathematical Problems in Engineering, 2015, 2015: 275057

    [16]

    Liu WC, Yao J, Chen ZX. Analytical solution of a double moving boundary problem for nonlinear flow in one-dimensional semi-infinite long porous media with low permeability. Acta Mechanica Sinica, 2014, 30(1): 50-58 doi: 10.1007/s10409-013-0091-5

    [17]

    Liu WC, Yao J, Wang YY. Exact analytical solutions of moving boundary problems of one-dimensional flow in semi-infinite long porous media with threshold pressure gradient. International Journal of Heat and Mass Transfer, 2012, 55(21-22): 6017-6022 doi: 10.1016/j.ijheatmasstransfer.2012.06.012

    [18] 朱维耀, 刘今子, 宋洪庆等.低/特低渗透油藏非达西渗流有效动用计算方法.石油学报, 2010, 31(3): 452-457 doi: 10.7623/syxb201003018

    Zhu Weiyao, Liu Jinzi, Song Hongqing, et al. Calculation of effective startup degree of non-Darcy flow in low or ultra-low permeability reservoirs. Acta Petrolei Sinica, 2010, 31(3): 452-457 (in Chinese) doi: 10.7623/syxb201003018

    [19]

    Zhu WY, Song HQ, Huang XH, et al. Pressure characteristics and effective deployment in a water-bearing tight gas reservoir with low-velocity non-Darcy flow. Energy & Fuels, 2011, 25(3): 1111-1117

    [20] 陈振宏, 王一兵, 郭凯等.高煤阶煤层气藏储层应力敏感性研究.地质学报, 2008, 82(10): 1390-1395. doi: 10.3321/j.issn:0001-5717.2008.10.013

    Chen Zhenhong, Wang Yibing, Guo Kai, et al. Stress sensitivity of high-rank coalbed methane reservoir. Acta Geologica Sinica, 2008, 82(10): 1390-1395 (in Chinese) doi: 10.3321/j.issn:0001-5717.2008.10.013

    [21]

    Liu WC, Liu YW, Niu CC, et al. Numerical investigations of a coupled moving boundary model of radial flow in low-permeable stress -sensitive reservoirs with threshold pressure gradient. Chinese Physics B, 2016, 25(2): 024701 doi: 10.1088/1674-1056/25/2/024701

    [22] 刘文超, 刘曰武, 门相勇等. 低渗煤层动边界不定常渗流模型建立//沈清主编. 第八届全国流体力学学术会议论文集. 北京: 科学出版社, 2014

    Liu Wenchao, Liu Yuewu, Men Xiangyong, et al. Construction of a model of unsteady seepage flow in low-permeable coalbed gas reservoirs with moving boundaries//Shen Qing, ed. Proceedings of the 8th National Conference on Fluid Mechanics. Beijing: Science Press, 2014 (in Chinese)

    [23]

    Liu WC, Liu YW, Niu CC, et al. A model of unsteady seepage flow in low-permeable coalbed with moving boundary in consideration of wellbore storage and skin effect. Procedia Engineering, 2015, 126: 517-521 doi: 10.1016/j.proeng.2015.11.294

    [24]

    Wan YZ, Liu YW, Ouyang WP, et al. Desorption area and pressure-drop region of wells in a homogeneous coalbed. Journal of Natural Gas Science and Engineering, 2016, 28: 1-14 doi: 10.1016/j.jngse.2015.11.026

    [25]

    Liu YW, Ouyang WP, Zhao PH, et al. Numerical well test for well with finite conductivity vertical fracture in coalbed. Applied Mathematics and Mechanics, 2014, 35(6): 729-740 doi: 10.1007/s10483-014-1825-6

    [26] 牛丛丛, 刘曰武, 蔡强等.煤层气井气水两相分布不稳定试井模型.力学与实践, 2013, 35(5): 35-41 doi: 10.6052/1000-0879-13-161

    Niu Congcong, Liu Yuewu, Cai Qiang, et al. Transient well test model for the well with gas and water distributed in coalbed. Mechanics in Engineering, 2013, 35(5): 35-41 (in Chinese) doi: 10.6052/1000-0879-13-161

    [27] 刘能强.实用现代试井解释方法.北京:石油工业出版社, 2008

    Liu Nengqiang. Practical Modern Well Testing Interpretation Method. Beijing: Petroleum Industry Press, 2008 (in Chinese)

    [28]

    Liu WC, Yao J, Chen ZX, et al. Effect of quadratic pressure gradient term on a one-dimensional moving boundary problem based on modified Darcy's law. Acta Mechanica Sinica, 2016, 32(1): 38-53 doi: 10.1007/s10409-015-0526-2

    [29] 刘文超, 姚军, 陈掌星等.低渗透多孔介质渗流动边界模型的解析与数值解.力学学报, 2015, 47(4): 605-612 doi: 10.6052/0459-1879-14-385

    Liu Wenchao, Yao Jun, Chen Zhangxin, et al. Research on analytical and numerical solutions of a moving boundary model of seepage flow in low -permeable porous media. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(4): 605-612 (in Chinese) doi: 10.6052/0459-1879-14-385

    [30] 李凡华, 刘慈群.含启动压力梯度的不定常渗流的压力动态分析.油气井测试, 1997, 6(1): 1-4 http://www.cnki.com.cn/Article/CJFDTOTAL-YQJC199701000.htm

    Li Fanhua, Liu Ciqun. Pressure transient analysis for unsteady porous flow with start-up pressure derivative. Well Testing, 1997, 6(1): 1-4 (in Chinese) http://www.cnki.com.cn/Article/CJFDTOTAL-YQJC199701000.htm

    [31]

    Clarkson CR, Behmanesh H, Chorney L. Production data and pressure transient analysis of Horsehoe Canyon CBM wells, Part Ⅱ: accounting for dynamic skin. SPE 148994, 2011

  • 期刊类型引用(8)

    1. 张鹏,王相春,封从军,郑力会,张妍,孙萌思. 基于多因素的煤储层气水两相非稳态流入动态评价方法及应用. 天然气地球科学. 2023(09): 1641-1651 . 百度学术
    2. 李东奇,杨志兵,张乐,胡冉,陈益峰. 空气-悬浮液驱替条件下颗粒边壁滞留研究. 力学学报. 2023(11): 2531-2538 . 本站查看
    3. 曹冲,程林松,张向阳,贾品,时俊杰. 基于多变量小样本的渗流代理模型及产量预测方法. 力学学报. 2021(08): 2345-2354 . 本站查看
    4. 姚成宝,付梅艳,韩峰,闫凯. 欧拉坐标系下具有锐利相界面的可压缩多介质流动数值方法研究. 力学学报. 2020(04): 1063-1079 . 本站查看
    5. 程志林,宁正福,曾彦,王庆,隋微波,张文通,叶洪涛,陈志礼. 格子Boltzmann方法模拟多孔介质惯性流的边界条件改进. 力学学报. 2019(01): 124-134 . 本站查看
    6. 张嫚嫚,孙姣,陈文义. 一种基于几何重构的Youngs-VOF耦合水平集追踪方法. 力学学报. 2019(03): 775-786 . 本站查看
    7. 王昭,严红. 基于气液相界面捕捉的统一气体动理学格式. 力学学报. 2018(04): 711-721 . 本站查看
    8. 牛丛丛,鞠斌山. 双重介质煤层气-水两相不稳定试井模型. 中国科技论文. 2018(03): 345-351+370 . 百度学术

    其他类型引用(5)

图(9)
计量
  • 文章访问数:  1445
  • HTML全文浏览量:  322
  • PDF下载量:  632
  • 被引次数: 13
出版历程
  • 收稿日期:  2017-01-09
  • 网络出版日期:  2017-04-09
  • 刊出日期:  2017-07-17

目录

/

返回文章
返回