裂缝型碳酸盐岩岩石物理模型的构建流程
(1)利用Voigt-Reuss-Hill或者Hashin-Shtrikman平均计算混合矿物的弹性模量;
(2)利用DEM模型将孔隙加进系统,并计算干燥岩石骨架的体积模量和剪切模量;
(3)利用Hudson理论和Schoenberg&Sayers理论,在碳酸盐岩介质各向同性背景中加入裂缝系统,修正计算出的模量;
(4)利用Wood公式将孔隙流体进行混合,计算出混合流体的体积模量;
( 5 )计算裂缝型碳酸盐岩的弹性参数和各向异性参数.
1.利用Voigt-Reuss-Hill或者Hashin-Shtrikman平均计算混合矿物的弹性模量
利用岩石基质组分,计算混合矿物的弹性模量
井的矿物成分含量和所需的井曲线如下:
![](https://img2023.cnblogs.com/blog/1602810/202301/1602810-20230112153449760-241269071.png)
![](https://img2023.cnblogs.com/blog/1602810/202301/1602810-20230112153720789-1244589856.png)
经过Hashin-Shtrikman计算后获得的等效矿物模型结果如下:
![](https://img2023.cnblogs.com/blog/1602810/202301/1602810-20230112153814249-482588535.png)
//////////////////////////////////////////////////////////////////////////////////////////////
void HashinShtrikman(float f[], int m, int n, float K[], float G[], float y[])//对应第一步 等效岩石基质的建立对应公式1和公式2
{
float *K_HS1 = new float[m];
float *G_HS1 = new float[m];
float *K_HS2 = new float[m];
float *G_HS2 = new float[m];
float Kmax = max(K, n);
float Kmin = min(K, n);
float Gmax = max(G, n);
float Gmin = min(G, n);
if (n == 2)
{
float K1 = K[0];
float K2 = K[1];
float G1 = G[0];
float G2 = G[1];
for (int i = 0; i < m; i++)
{
float f1 = f[2 * i];
float f2 = f[2 * i + 1];//组分体积比归一化
K_HS1[i] = K1 + f2 / (1.0f / (K2 - K1) + f1 / (K1 + 4.0f / 3.0f * Gmin));
G_HS1[i] = G1 + f2 / (1.0f / (G2 - G1) + f1 / (G1 + Gmin / 6 * (9 * Kmin + 8 * Gmin) / (Kmin + 2 * Gmin))); //公式1
K_HS2[i] = K1 + f2 / (1.0f / (K2 - K1) + f1 / (K1 + 4.0f / 3.0f * Gmax));
G_HS2[i] = G1 + f2 / (1.0f / (G2 - G1) + f1 / (G1 + Gmax / 6 * (9 * Kmax + 8 * Gmax) / (Kmax + 2 * Gmax))); //公式2
y[4 * i + 0] = K_HS1[i];
y[4 * i + 1] = K_HS2[i];
y[4 * i + 2] = G_HS1[i];
y[4 * i + 3] = G_HS2[i];
}
}
else
{
//float Kmax = max(K, n);
//float Kmin = min(K, n);
//float Gmax = max(G, n);
//float Gmin = min(G, n);
float Emax = Gmax / 6.0f * (9 * Kmax + 8 * Gmax) / (Kmax + 2 * Gmax);
float Emin = Gmin / 6.0f * (9 * Kmin + 8 * Gmin) / (Kmin + 2 * Gmin);
for (int i = 0; i < m; i++)
{
K_HS1[i] = 0;
K_HS2[i] = 0;
G_HS1[i] = 0;
G_HS2[i] = 0;
for (int j = 0; j < n; j++)//
{
K_HS1[i] = K_HS1[i] + f[i * n + j] / (K[j] + 4.0f / 3.0f * Gmax);//
K_HS2[i] = K_HS2[i] + f[i * n + j] / (K[j] + 4.0f / 3.0f * Gmin);
G_HS1[i] = G_HS1[i] + f[i * n + j] / (G[j] + Emax);
G_HS2[i] = G_HS2[i] + f[i * n + j] / (G[j] + Emin);
}
K_HS1[i] = 1.0f / K_HS1[i] - 4.0f / 3.0f * Gmax;
K_HS2[i] = 1.0f / K_HS2[i] - 4.0f / 3.0f * Gmin;
G_HS1[i] = 1.0f / G_HS1[i] - Emax;
G_HS2[i] = 1.0f / G_HS2[i] - Emin;
y[4 * i + 0] = K_HS1[i];
y[4 * i + 1] = K_HS2[i];
y[4 * i + 2] = G_HS1[i];
y[4 * i + 3] = G_HS2[i];
}
}
delete[]K_HS1;
delete[]K_HS2;
delete[]G_HS1;
delete[]G_HS2;
}
2.利用DEM模型将孔隙加进系统,并计算干燥岩石骨架的体积模量和剪切模量;
孔隙空间分成主成分孔隙和次成分孔隙
![](https://img2023.cnblogs.com/blog/1602810/202301/1602810-20230112154437741-1938587576.png)
![](https://img2023.cnblogs.com/blog/1602810/202301/1602810-20230112155151842-1410336431.png)
利用广义xuwhite模型计算岩石的弹性模量:
![](https://img2023.cnblogs.com/blog/1602810/202301/1602810-20230112155228359-1438232232.png)
![](https://img2023.cnblogs.com/blog/1602810/202301/1602810-20230112155304279-1626186451.png)
![](https://img2023.cnblogs.com/blog/1602810/202301/1602810-20230112155331852-1379575340.png)
![](https://img2023.cnblogs.com/blog/1602810/202301/1602810-20230112155356080-273002535.png)
其中:
![](https://img2023.cnblogs.com/blog/1602810/202301/1602810-20230112155429382-864358024.png)
采用微分等效介质模型加入孔隙。
初始条件是
K(0)=K1和μ(0)=μ1;
其中:
K1,μ1为初始主相材料的体积模量和剪切模量
K2,μ2为逐渐加入包含物材料的体积模量和剪切模量
y=相2的含量
对于流体包含物和空包含物,y等于孔隙度φ,项P和Q是给定的几何因素,P和Q的上标*2指的是几何因素具有等效的背景介质中的包含物材料(具体背景材料参数和系数参考岩石物理手册P110)。
DEM计算后的结果:
![](https://img2023.cnblogs.com/blog/1602810/202301/1602810-20230112155459982-650888970.png)
void XuWhite_DEM2(float Por, float Sg, float So, float Vsh, int MatrixMixMethod, int PorGeometry, int NAspect, float *y)//DEM微分等效模型公式3
{
float Km = 71.07;
float Gm = 28.3833;
float Pm = 2.71;
float Ksh = 42.05;
float Gsh = 19.08;
float Psh = 2.62;
float Kw = 2.25;
float Gw = 0;
float Pw = 1;
float Kg = 0.71415;
float Gg = 0;
float Pg = 0.54;
float Ko = 0.91415;
float Po = 0.77;
float ratio[2] = { 0.12,0.05 };
float Vsand, Vsand1, Vsh1, Por_s, Por_c;
//matrix
Vsand = 1 - Vsh - Por;
Vsh1 = Vsh / (1 - Por);
Vsand1 = Vsand / (1 - Por);
Por_s = Vsand1 * Por;//体积比 Por_s 对应Tiijj0[0],Por_c对应Tiijj[1]即:组分1为砂岩,组分2为泥岩
Por_c = Vsh1 * Por;
//float f[2] = { 1 - Vsh, Vsh };//这个参数传递给HashinShtrikman()
float f[2] = { Vsand1, Vsh1 };
float K[2] = { Km, Ksh };
float G[2] = { Gm, Gsh };
float P[2] = { Pm,Psh };
float *y1 = new float[4];
////*********Dry matrix model (using DEM model)
// Xu-white coefficient
//float A,B,R,Tiijj,FF,alpha,theta,f,F1,F2,F3,F4,F5,F6,F7,F8,F9;
int NP = 5000;//NP是积分的个数,确定积分间隔
// int N=length(Fraction);
float Fraction[2] = { Por_s, Por_c };//第一孔隙是主岩性,第二孔隙是副岩性:例如砂泥组合第一孔隙是砂岩,第二孔隙是泥岩
float Aspect[2] = { ratio[0],ratio[1] };//第一孔隙是主岩性,第二孔隙是副岩性;
float p, dp;
float PP, QQ, Em, bm;
float K_in, Mu_in, Ein;
float K_m, Mu_m, K_tmp, Mu_tmp, Kd, Gd;
int i, j;
K_in = 0;//填充空气
Mu_in = 0;
K_m = Km;
Mu_m = Gm;
p = 0;
//for(i=0;i<NAspect;i++)//几种孔隙
for (i = NAspect - 1; i >= 0; i--)//几种孔隙
{
dp = Fraction[i] / NP;
K_tmp = K_m;
Mu_tmp = Mu_m;
for (j = 1; j <= NP; j++)
{
float A = -1;
float B = 0;
//float R = 3 * Gm / (3 * Km + 4 * Gm);
//float *Tiijj= new float[2];
//float *FF = new float[2];
float R = 3 * Mu_tmp / (3 * K_tmp + 4 * Mu_tmp);//K_tmp,Mu_tmp是每次的增加孔隙中间结果
float Tiijj, FF;//不需要记录每个孔隙
//for (int i = 0; i < 2; i++)
//{
float alpha = Aspect[i];
float theta = alpha / powf(1 - alpha * alpha, 1.5) * (acos(alpha) - alpha * sqrt(1 - alpha * alpha));
float f = alpha * alpha / (1 - alpha * alpha) * (3 * theta - 2);
float F1 = 1 + A * (1.5 * (f + theta) - R * (1.5 * f + 2.5 * theta - 4.0 / 3.0));
float F2 = 1 + A * (1 + 1.5 * (f + theta) - R / 2 * (3 * f + 5 * theta)) + B * (3 - 4 * R) + A / 2 * (A + 3 * B) * (3 - 4 * R) * (f + theta - R * (f - theta + 2 * theta * theta));
float F3 = 1 + A * (1 - (f + 1.5 * theta) + R * (f + theta));
//float F3 = 1 +A/2*( R*(2-theta)+(1+alpha*alpha)/alpha*f*(R-1) );
float F4 = 1 + A / 4 * (f + 3 * theta - R * (f - theta));
float F5 = A * (-f + R * (f + theta - 4.0 / 3.0)) + B * theta * (3 - 4 * R);
float F6 = 1 + A * (1 + f - R * (f + theta)) + B * (1 - theta) * (3 - 4 * R);
float F7 = 2 + A / 4 * (3 * f + 9 * theta - R * (3 * f + 5 * theta)) + B * theta * (3 - 4 * R);
float F8 = A * (1 - 2 * R + f / 2 * (R - 1) + theta / 2 * (5 * R - 3)) + B * (1 - theta) * (3 - 4 * R);
float F9 = A * ((R - 1) * f - R * theta) + B * theta * (3 - 4 * R);
Tiijj = 3 * F1 / F2;
FF = 2 / F3 + 1 / F4 + (F4 * F5 + F6 * F7 - F8 * F9) / (F2 * F4);
//printf("======%d=====%d=====%f=====%f===\n",i,j,alpha,theta,f);
//}
PP = 1.0 / 3.0*Tiijj;
QQ = 1.0 / 5.0*FF;
p = p + dp;
//p=dp;
K_tmp = (K_in - K_tmp)*PP*dp / (1 - p) + K_tmp;
Mu_tmp = (Mu_in - Mu_tmp)*QQ*dp / (1 - p) + Mu_tmp;
}
K_m = K_tmp;
Mu_m = Mu_tmp;
}
Kd = K_m;
Gd = Mu_m;
y[0] = Kd;
y[1] = Gd;
y[2] = Pm;
}
3.利用Hudson理论和Schoenberg & Sayers理论,在碳酸盐岩介质各向同性背景中加入裂缝系统,修正计算出的模量。
从前面1,2步建立的等效岩石背景基质中,获得基质的体积模量和剪切模量。在其基础上加入裂缝的影响,修正计算的模量。
目前有两种应用广泛的裂缝模型,Hudson 模型和Schoenberg线性滑动模型。Hudson理论是建立在整体平均的基础上。该模型假设弹性介质的内部分布着薄硬币形状的椭球缝隙。Hudson模型利用裂缝密度和缝隙纵横比来描述裂缝系统的结构。裂缝密度e的定义如下:
其中,φ代表孔隙度,α 是裂隙纵横比.则Hudson模型的弹性矩阵犆C表示为:
式中,C0为各向同性背景弹性矩阵,C1,C2为引入裂隙的各向异性一阶校正和二阶校正。
![](https://img2023.cnblogs.com/blog/1602810/202301/1602810-20230112155557052-344713635.png)
![](https://img2023.cnblogs.com/blog/1602810/202301/1602810-20230112155628784-779493318.png)
Schoenberg模型的弹性矩阵由各向同性参数λ 和μ,以及垂直于裂缝面和平行于裂缝面的弹性差值ΔN 和ΔT 进行描述。表示裂缝对骨架岩石各向同性特征的影响。在裂缝诱导等效各向异性介质中,裂缝等效介质弹性参数矩阵C****表示为背景各向同性介质弹性参数Cb与裂缝介质弹性参数Cf之差;Schoenberg线性滑动模型的弹性矩阵为:
![](https://img2023.cnblogs.com/blog/1602810/202301/1602810-20230112155655818-1432944745.png)
△N和△T分别表示裂缝的法向弱度和切向弱度
下面列出Schoenberg模型的计算代码:
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////计算裂缝密度,ΔN 和ΔT
for (int i = 0; i < nsample; i++)
{
Sg[i] = 1 - Sw[i];
e[i] = 3 * FraPOR[i] / (4 * PI*arf);
}
for (int i = 0; i < nsample; i++)
{
g[i] = (Vs_model[i] * Vs_model[i] / Vp_model[i] / Vp_model[i]);
deltaN[i] = 4 * e[i] / (3 * g[i] * (1 - g[i])*(1 + (k_f + 1.3333*miu_f) / (PI*(1 - g[i])*miu[i] * arf)));
deltaT[i] = 16 * e[i] / (3 * (3 - 2 * g[i])*(1 + 4 * miu_f / (PI*(3 - 2 * g[i])*miu[i] * arf)));
}
4.利用Wood公式将孔隙流体进行混合,计算出混合流体的体积模量;
wood公式进行的平均等效模型表达式为:
其中f,K,ρ分别为成分的体积含量,体积模量和密度;
// wood流体混合模式
{
if (nFluidCom == 0)
{
Kf = 1 / (So / Ko + (1 - So) / Kw);//水油
}
if (nFluidCom == 1)
{
Kf = 1 / (Sg / Kg + (1 - Sg) / Kw);//水气
}
else
{
Kf = 1 / (So / Ko + Sg / Kg + (1 - Sg - So) / Kw);//气油水
}
}
( 5 )计算裂缝型碳酸盐岩的弹性参数和各向异性参数.
由于研究区域只含有气,不用进行流体替换。
各向异性参数的表达式为:
对应的计算代码:
epsilon[i] = (CHTI[0][0] - CHTI[2][2]) / 2 / CHTI[2][2];
delta[i] = ((CHTI[0][2] + CHTI[4][4])*(CHTI[0][2] + CHTI[4][4]) - (CHTI[2][2] - CHTI[4][4])*(CHTI[2][2] - CHTI[4][4])) / (2 * CHTI[2][2] * (CHTI[2][2] - CHTI[4][4]));
gamma[i] = (CHTI[5][5] - CHTI[3][3]) / (2 * CHTI[3][3]);
![](https://img2023.cnblogs.com/blog/1602810/202301/1602810-20230112155742202-77471972.png)
计算所得的各向异性参数与实测的吻合度较高。
参考及其引用文献:
张广智,陈怀震,王琪等.基于碳酸盐岩裂缝岩石物理模型的横波速度和各向异性参数预测.地球物理学报,2013,56(5):1707
1715,doi:10.6038/cjg20130528.
岩石物理手册.Gray Mavko and so on.