CN106153859A - 一种盐岩动水溶蚀的试验装置和计算方法 - Google Patents

Abstract

本发明公开了一种盐岩动水溶蚀的试验装置和计算方法。该试验装置简单，所采用的试验组成设备为常见装置，通过两个接头将通水软管和水龙头连接，所用试验组成设备无特殊要求。根据盐岩溶解特性和渗透性的耦合因素，建立了盐岩动水溶蚀模型的常微分方程组，并给出了该动水溶蚀模型的计算方法，通过建立有限差分网格，运用有限差分法数值求解盐岩动水溶蚀模型的数值解，引入粒子群算法对所建立的模型进行了参数反演，反演可以得到不同流量条件下的盐岩动水溶蚀参数。解决了在盐岩动水溶蚀试验等方面的试验难题和理论计算难题，具有良好的理论价值和实际意义。

Description

一种盐岩动水溶蚀的试验装置和计算方法

技术领域

本发明涉及矿业工程领域中盐岩水溶开采，具体涉及一种盐岩动水溶蚀的试验装置和计算方法。

背景技术

水溶开采技术始于盐岩采卤，现已成为一门新兴的开采方法。它不仅广泛应用于盐岩水溶采卤，还应用于一些金属矿物的开采。近年来，由于盐岩矿床采卤后溶腔的利用，在地下深处建造大量的盐腔，存储碳氢化合物和处置化学工业的有害物质，以及核工业的放射性废料，使之成为专门建造地下盐腔储库的新方法和新技术，并已发展成为除盐岩生产外，包括石油、能源、化工、核工业、环境保护以及国防等多部门所共同关注和研究的技术。

受试验条件的限制，以往关于盐岩动水溶蚀特性方面的试验装置很少见。有鉴于此，需要发明一种盐岩水溶蚀的试验装置和计算方法，用以为盐岩动水溶蚀机理的研究提供试验方法和计算依据。

发明内容

本发明的目的在于克服现有技术的缺陷，提供一种岩动水溶蚀特性的试验装置，并建立了相应的盐岩的动水溶蚀模型，给出了动水溶蚀模型的计算方法。

为了实现上述目的，本发明一方面采取的技术方案如下：

一种盐岩动水溶蚀的试验装置，其包括圆柱形盐岩试件1、用于夹持固定盐岩试件1的夹持装置和溶蚀水供应装置；

所述盐岩试件轴向开设有通水小孔；

所述溶蚀水供应装置通过软管和所述盐岩试件的通水小孔相连通。

进一步地，所述夹持装置包括铁架台和固定于所述铁架台上的烧杯链夹，所述盐岩试件通过烧杯链夹进行固定。

进一步地，所述软管上安装有液体流量计。

进一步地，所述盐岩试件的外部涂覆有防水涂层。

本发明另一方面提供了一种盐岩动水溶蚀的计算方法，具体如下：

依据上述实验装置测量不同流量条件下的不同时间对应的盐岩试件1的质量损失，得到不同流量条件下盐岩的动水溶蚀模型的常微分方程组，运用有限差分法数值求解盐岩动水溶蚀模型的数值解，引入粒子群算法对所建立的模型进行了参数反演，反演得到不同流量条件下的盐岩动水溶蚀参数。

进一步地，所述的动水溶蚀模型的常微分方程组如下：

$\left\{\begin{array}{c}\frac{\∂R}{\∂t}=\frac{2}{\mathrm{\δ}}D\frac{M}{{\mathrm{\ρ}}_s}(C_s-C_{xt})\\ \frac{\∂C}{\∂x}=\frac{2D}{\mathrm{\δ}Q}\mathrm{\π}(d+2R_{xt})(C_s-C_{xt})\\ R_{xt}{|}_{t=0}=0\\ C_{xt}{|}_{x=0}=C_0\end{array}\right.---\left(10\right)$

式中：

ρ_s为盐岩的密度，单位为cm³·s^-1；

M为盐岩溶液的摩尔质量，单位为g·mol^-1；

δ为边界层厚度，单位为cm；

C_s为盐岩溶液的饱和浓度，单位为mol·L-¹；

C_xt为t时刻在孔隙x处的溶液浓度，单位为mol·L^-1；

R为盐岩固壁的溶解厚度，单位为cm；

C为边界层内溶液浓度，单位为mol·L^-1；

D为扩散系数，单位为cm²·s^-1；

R_xt为在t时刻x位置处的溶解厚度，单位为cm；

d为小孔初始直径，单位为cm；

Q为流量大小，单位为cm³·s^-1；

C₀为溶液初始浓度，单位为mol·L^-1。

进一步地，所述采用粒子群算法对所建立的模型的参数反演中，所建立的盐岩质量与溶蚀半径的计算公式：

$R\left(k\right)=\sqrt{\frac{m\left(k\right)}{{\mathrm{\ρ}}_{s}l\mathrm{\π}}+\frac{d^2}{4}}-\frac{d}{2}---\left(15\right)$

式中：

R(k)为k时刻盐岩溶蚀半径，单位为cm；

l为试样高度，单位为cm；

m(k)为t时刻溶蚀质量，单位为g；

d为小孔初始直径，单位为cm；

ρ_s为盐岩的密度，单位为cm³·s^-1。

进一步地，所述采用粒子群算法对所建立的模型的参数反演中，所建立的适应度函数如下：

$f=\frac{1}{N}\underset{k=1}{\overset{N}{\Σ}}(\sqrt{\frac{m\left(k\right)}{{\mathrm{\ρ}}_{s}l\mathrm{\π}}+\frac{d^2}{4}}-\frac{d}{2}-R)---\left(16\right)$

式中：

R(k)为k时刻盐岩溶蚀半径，单位为cm；

l为试样高度，单位为cm；

m(k)为t时刻溶蚀质量，单位为g；

d为小孔初始直径，单位为cm；

ρ_s为盐岩的密度，单位为cm³·s^-1；

R为盐岩固壁的溶解厚度，单位为cm；

N为试验记录次数。

进一步地，所述用粒子群算法对所建立的模型的参数反演中，基于粒子群算法编写的主程序，命名为main.m，该程序主要步骤包括：

调用cc.Mat文件，获取随机位移矩阵和速度矩阵，然后执行Numerical.m文件，计算出溶蚀半径R，并提取计算出的结果带入适应度函数，得到每个粒子的适应度大小；

2)对每个粒子适应度值进行比较，找到全局最优解(即最小值)，记录全局最优值的数值大小和位置；

3)按照粒子群算法原理更新速度矩阵和位置矩阵，并将更新的数据保存覆盖第一部分cc.Mat文件；

4)得到新的粒子速度和位置矩阵后，再次执行Numerical.m文件，计算出溶蚀半径R，并得到每个粒子的适应度值；

5)对速度、位置更新前后的每个粒子适应度值进行比较，如果更新后的适应度值更小，则保存记录该适应度值大小和位置，否则，则不保存；

6)比较全局最优值与速度、位置更新后的每个粒子的适应度值，如果更新后的某个粒子适应度值更小，则保存记录该适应度值大小和位置，否则，则不保存；

7)重复第4-6步，直到满足求解条件。

本发明中，所述公式(10)的推导原理和过程如下：

为了分析计算的方便，考虑理想条件下的盐岩试件，并作如下简化：①盐岩均质，溶解特性为各向同性，且无夹层；②动水溶解中的渗流为层流运动；③不可溶杂质随流体流动，在孔隙中不沉淀；④忽略不可溶杂质对盐溶液扩散的影响；⑤不考虑温差及温度变化的影响。

对于微元体dx，根据物质平衡原理并参照附图2可得:

$J\·d\·dx\·dt=\frac{{\mathrm{\ρ}}_s}{M}d\·dR\·dx---\left(1\right)$

式中：ρ_s为盐岩的密度，单位为cm³·s^-1；

M为盐岩溶液的摩尔质量，单位为g·mol^-1；

R为盐岩固壁的溶解厚度，单位为cm；

d为小孔直径，单位为cm。

盐溶边界层浓度剖面的分布可看作抛物线的形状，则溶蚀边界层浓度分布可用下式表示：

$C-C_1=(C_0-C_1){(1-\frac{z}{\mathrm{\δ}})}^2---\left(2\right)$

式中：

z为距盐岩壁面的距离，单位为cm；

C₀为盐岩壁面浓度，单位为cm；

C₁为边界层以外溶液平均浓度，单位为mol·L^-1；

C为边界层内溶液浓度，单位为mol·L^-1；

边界层厚度为δ，单位为cm。

可以建立盐岩固壁的溶解速度方程：

$\frac{dR}{dt}=-D\frac{M}{{\mathrm{\ρ}}_s}\frac{\∂C}{\∂n}{|}_{{\mathrm{\Γ}}_1}---\left(3\right)$

式中：Г₁为溶解接触面。

对公式(2)求导数，并假设扩散过程瞬间完成，整理后可得：

$\frac{\∂C}{\∂n}{|}_{{\mathrm{\Γ}}_1}=-\frac{2}{\mathrm{\δ}}(C_s-C_{xt})---\left(4\right)$

式中：C_s为盐岩溶液的饱和浓度(mol·L^-1)，C_xt为t时刻在孔隙x处的溶液浓度，其余符号意义同前(mol·L^-1)。

将公式(4)带入公式(3)可得：

$\frac{dR}{dt}=\frac{2}{\mathrm{\δ}}D\frac{M}{{\mathrm{\ρ}}_s}(C_s-C_{xt})---\left(5\right)$

盐岩溶解过程还受到流体的流速v和孔隙直径d控制，为此，需要建立盐岩溶解的平衡方程。设m_x为dt时间内流经x处孔隙断面的溶液离子摩尔数，m_x+dx为dt时间内流经x+dx处孔隙断面的溶液离子摩尔数，d为孔隙的初始直径，单位为cm，则有：

$\left\{\begin{array}{c}{m}_x=C_{xt}\frac{\mathrm{\π}}{4}{(d+2R_{xt})}^2{v}_{xt}dt\\ m_{x+dx}=(C_{xt}+\frac{\∂C}{\∂x}dx)\frac{\mathrm{\π}}{4}{(d+2R_{xt})}^2{v}_{xt}dt\end{array}\right.---\left(6\right)$

根据物质平衡原理，忽略盐岩溶解引起的孔隙内流体体积的变化，m_x+dx与m_x之差值即为dt内dx段溶解的盐岩摩尔数，即

$m_{x+dx}-m_x=\frac{{\mathrm{\ρ}}_s}{M}\mathrm{\π}(d+2R_{xt})\∂R\∂x---\left(7\right)$

假设不计流量损失，则有：

$Q=\frac{\mathrm{\π}}{4}{(d+2R_{xt})}^2{v}_{xt}---\left(8\right)$

式中：Q为流量大小，单位为cm³·s^-1。

公式(6)、(7和(8)联立，整理后可得：

$\frac{\∂C}{\∂x}=\frac{{\mathrm{\ρ}}_s}{MQ}\mathrm{\π}(d+2R_{xt})\frac{\∂R}{\∂t}---\left(9\right)$

由公式(5)和(9)构成偏微分方程组：

$\left\{\begin{array}{c}\frac{\∂R}{\∂t}=\frac{2}{\mathrm{\δ}}D\frac{M}{{\mathrm{\ρ}}_s}(C_s-C_{xt})\\ \frac{\∂C}{\∂x}=\frac{2D}{\mathrm{\δ}Q}\mathrm{\π}(d+2R_{xt})(C_s-C_{xt})\\ R_{xt}{|}_{t=0}=0\\ C_{xt}{|}_{x=0}=C_0\end{array}\right.---\left(10\right)$

与现有技术相比，本发明所取得的有益效果如下：

本发明所设计研制的盐岩动水溶蚀试验装置，结构简单，所采用的试验组成设备为常见装置，通过两个接头将通水软管和水龙头连接，所用试验组成设备无特殊要求。

可以对不同流量下的盐岩动水溶蚀过程展开试验研究，并建立了相应的盐岩的动水溶蚀模型，给出了有限差分法数值求解动水溶蚀模型的计算方法，本发明可为盐岩动水溶蚀机理的研究提供试验方法和计算依据，该试验装置和计算方法解决了在盐岩动水溶蚀试验等方面的试验难题和理论计算难题，具有良好的理论价值和实际意义。

附图说明

附图1为本发明的结构示意图；

附图2为本发明中盐岩动水溶解关系示意图

附图3为差分网格示意图；

附图4为粒子群算法反演程序流程图。

在附图中：

1盐岩试件、2软管、3铁架台、4烧杯链夹、5阀门、6液体流量计。

具体实施方式

以下结合附图对本发明进行进一步详细的叙述。

1、实验装置

如附图1所述的一种盐岩动水溶蚀的试验装置，其包括外部涂覆有防水涂层圆柱形盐岩试件1、用于夹持固定盐岩试件1的夹持装置和溶蚀水供应装置，所述夹持装置包括铁架台3和固定于所述铁架台3上的烧杯链夹4，所述盐岩试件1通过烧杯链夹4进行固定。

在盐岩试件1轴向开设有通水小孔，溶蚀水供应装置通过软管2和所述盐岩试件1的通水小孔相连通，所述软管2上安装有液体流量计6用于实时控制动水的流量，为了及时的控制开关，在软管2上设置有阀门5。

本实施例中盐岩试样为50×100mm的圆柱试件，通水小孔的直径为0.6cm。

2、试验方法如下：

①用电子天平对已经钻孔好的试件称重，用游标卡尺测量试件的长度和直径；

②用防水涂料703胶将试件密封，密封好后用电子天平称重；

③通过调整溶蚀水供应装置水龙头阀门开关来调节自来水流量使流量计示数为20L/h；

④将涂有防水涂料密封好的试件用烧杯链夹固定好，将通水软管插入试件中进行盐岩动水试验，并开始计时；

⑤每隔3min，拔出通水软管，并为保证流量的稳定将软管放入另一个大烧杯中，取下试件，擦干、吹风机吹干，并用电子天平称重；

⑥重复⑤步试验步骤直到溶解到试件的外表面；

⑦调节流量大小依次为30L/h、40L/h、50L/h，重复第④、⑤、⑥步试验步骤。

3、数值计算盐岩动水溶蚀模型

在不同流量条件下的盐岩动水溶蚀模型中，由于速度场和浓度场的相互影响，偏微分方程组(10)具有非线性特征，非常难以求出解析解，因此可以通过合适的离散格式，将偏微分方程组(10)离散化，通过数值方法进行求解计算。

可建立差分网格如图3，其中x的取值范围是[0,10](单位为cm)，t的取值范围是[0,2340](单位为s)。其中x的步长为h，t的步长为τ。取x＝nh,t＝iτ，并设

$\left\{\begin{array}{c}R(x,t)=R(nh,i\mathrm{\τ})=R_i^n\\ C(x,t)=C(nh,i\mathrm{\τ})=C_i^n\\ \frac{\∂R(x,t)}{\∂t}=\frac{\∂R}{\∂t}\\ \frac{\∂C(x,t)}{\∂x}=\frac{\∂C}{\∂x}\end{array}\right.$

由公式(10)建立差分格式如下：

$\left\{\begin{array}{c}\frac{R(nh,(i+1\left)\mathrm{\τ}\right)-R(nh,i\mathrm{\τ})}{h}=\frac{2}{\mathrm{\δ}}D\frac{M}{{\mathrm{\ρ}}_s}(C_s-C(nh,i\mathrm{\τ}\left)\right)\\ \frac{C\left(\right(n+1)h,i\mathrm{\τ})-C(nh,i\mathrm{\τ})}{h}=\frac{2D}{\mathrm{\δ}Q}\mathrm{\π}(d+2R(nh,i\mathrm{\τ}\left)\right)(C_s-C(nh,i\mathrm{\τ}\left)\right)\end{array}\right.---\left(11\right)$

令可得：

$\left\{\begin{array}{c}{R}_{i+1}^n=\frac{D}{\mathrm{\δ}}(a-{\mathrm{bC}}_i^n)+R_i^n\\ C_i^{n+1}=\frac{D}{\mathrm{\δ}}(A-{\mathrm{BC}}_i^n)(d+2R_i^n)+C_i^n\end{array}\right.---\left(12\right)$

可得：

$\left\{\begin{array}{c}{R}_i^n=\frac{D}{\mathrm{\δ}}(a-{\mathrm{bC}}_{i-1}^n)+R_{i-1}^n\\ C_i^n=\frac{D}{\mathrm{\δ}}(A-{\mathrm{BC}}_i^{n-1})(d+2R_i^{n-1})+C_i^{n-1}\end{array}\right.---\left(13\right)$

进一步可得：

$\left\{\begin{array}{c}{R}_i^0=\frac{D}{\mathrm{\δ}}(a-{\mathrm{bC}}_{i-1}^0)+R_{i-1}^0\\ C_0^n=\frac{D}{\mathrm{\δ}}(A-{\mathrm{BC}}_0^{n-1})(d+2R_0^{n-1})+C_0^{n-1}\end{array}\right.---\left(14\right)$

由初始条件R(x,0)＝0,C(0,t)＝C₀可知由于是已知的，且由递推公式(14)可得进而初始条件为已知。

4、粒子群算法参数反演

反演程序在MATLAB软件中实现，程序流程图如图4所示。一共包括三个部分，第一部分是速度矩阵和位置矩阵的随机初始化Mat文件，命名为cc.Mat；

第二部分是运用有限差分法求解溶蚀半径的M程序，命名为Numerical.m，该程序主要是调用第一部分生成的随机速度矩阵和随机位移矩阵来计算溶蚀半径R；

第三部分是基于粒子群算法编写的主程序，命名为main.m，该程序主要步骤包括：

调用cc.Mat文件，获取随机位移矩阵和速度矩阵，然后执行Numerical.m文件，计算出溶蚀半径R，并提取计算出的结果带入适应度函数，得到每个粒子的适应度大小；

2)对每个粒子适应度值进行比较，找到全局最优解(即最小值)，记录全局最优值的数值大小和位置；

3)按照粒子群算法原理更新速度矩阵和位置矩阵，并将更新的数据保存覆盖第一部分cc.Mat文件；

4)得到新的粒子速度和位置矩阵后，再次执行Numerical.m文件，计算出溶蚀半径R，并得到每个粒子的适应度值；

5)对速度、位置更新前后的每个粒子适应度值进行比较，如果更新后的适应度值更小，则保存记录该适应度值大小和位置，否则，则不保存；

6)比较全局最优值与速度、位置更新后的每个粒子的适应度值，如果更新后的某个粒子适应度值更小，则保存记录该适应度值大小和位置，否则，则不保存；

7)重复第4-6步，直到满足求解条件。

以上所述实施方式仅为本发明的优选实施例，而并非本发明可行实施的穷举。对于本领域一般技术人员而言，在不背离本发明原理和精神的前提下对其所作出的任何显而易见的改动，都应当被认为包含在本发明的权利要求保护范围之内。

Claims (8)

1.一种盐岩动水溶蚀的试验装置，其特征在于，其包括圆柱形盐岩试件(1)、用于夹持固定盐岩试件(1)的夹持装置和溶蚀水供应装置；

所述盐岩试件(1)轴向开设有通水小孔；

所述溶蚀水供应装置通过软管(2)和所述盐岩试件(1)的通水小孔相连通。

2.根据权利要求1所述的一种盐岩动水溶蚀的试验装置，其特征在于，所述夹持装置包括铁架台(3)和固定于所述铁架台(3)上的烧杯链夹(4)，所述盐岩试件(1)通过烧杯链夹(4)进行固定。

3.根据权利要求1所述的一种盐岩动水溶蚀的试验装置，其特征在于，所述软管(2)上安装有液体流量计(6)。

4.根据权利要求1所述的一种盐岩动水溶蚀的试验装置，其特征在于，所述盐岩试件(1)的外部涂覆有防水涂层。

5.一种盐岩动水溶蚀的计算方法，其特征在于，具体如下：

依据权利1-4任一项所述的实验装置测量不同流量条件下的不同时间对应的盐岩试件的质量损失，得到不同流量条件下盐岩的动水溶蚀模型的常微分方程组，运用有限差分法数值求解盐岩动水溶蚀模型的数值解，引入粒子群算法对所建立的模型进行了参数反演，反演得到不同流量条件下的盐岩动水溶蚀参数。

6.根据权利要求5所示的一种盐岩动水溶蚀的计算方法，其特征在于，所述的动水溶蚀模型的常微分方程组如下：

$\left\{\begin{array}{c}\frac{\∂R}{\∂t}=\frac{2}{\mathrm{\δ}}D\frac{M}{{\mathrm{\ρ}}_s}(C_s-C_{xt})\\ \frac{\∂C}{\∂x}=\frac{2D}{\mathrm{\δ}Q}\mathrm{\π}(d+2R_{xt})(C_s-C_{xt})\\ R_{xt}{|}_{t=0}=0\\ C_{xt}{|}_{x=0}=C_0\end{array}\right.---\left(10\right)$

式中：

ρ_s为盐岩的密度，单位为cm³·s^-1；

M为盐岩溶液的摩尔质量，单位为g·mol^-1；

δ为边界层厚度，单位为cm；

C_s为盐岩溶液的饱和浓度，单位为mol·L^-1；

C_xt为t时刻在孔隙x处的溶液浓度，单位为mol·L^-1；

R为盐岩固壁的溶解厚度，单位为cm；

C为边界层内溶液浓度，单位为mol·L^-1；

D为扩散系数，单位为cm²·s^-1；

R_xt为在t时刻x位置处的溶解厚度，单位为cm；

d为小孔初始直径，单位为cm；

Q为流量大小，单位为cm³·s^-1；

C₀为溶液初始浓度，单位为mol·L^-1。

7.根据权利要求5所示的一种盐岩动水溶蚀的计算方法，其特征在于，所述采用粒子群算法对所建立的模型的参数反演中，所建立的盐岩质量与溶蚀半径的计算公式：

$R\left(k\right)=\sqrt{\frac{m\left(k\right)}{{\mathrm{\ρ}}_{s}l\mathrm{\π}}+\frac{d^2}{4}}-\frac{d}{2}---\left(15\right)$

式中：

R(k)为k时刻盐岩溶蚀半径，单位为cm；

l为试样高度，单位为cm；

m(k)为t时刻溶蚀质量，单位为g；

d为小孔初始直径，单位为cm；

ρ_s为盐岩的密度，单位为cm³·s^-1。

8.根据权利要求5所示的一种盐岩动水溶蚀的计算方法，其特征在于，所述采用粒子群算法对所建立的模型的参数反演中，所建立的适应度函数如下：

$f=\frac{1}{N}\underset{k=1}{\overset{N}{\Σ}}(\sqrt{\frac{m\left(k\right)}{{\mathrm{\ρ}}_{s}l\mathrm{\π}}+\frac{d^2}{4}}-\frac{d}{2}-R)---\left(16\right)$

式中：

R(k)为k时刻盐岩溶蚀半径，单位为cm；

l为试样高度，单位为cm；

m(k)为t时刻溶蚀质量，单位为g；

d为小孔初始直径，单位为cm；

ρ_s为盐岩的密度，单位为cm³·s^-1；

R为盐岩固壁的溶解厚度，单位为cm；

N为试验记录次数。