首页 > 其他分享 >裂缝型碳酸盐岩岩石物理模型的构建流程

裂缝型碳酸盐岩岩石物理模型的构建流程

时间：2023-01-12 16:13:12浏览次数：52  
标签：CHTI 流程 float Por HS2 HS1 碳酸盐岩 裂缝

（１）利用Voigt-Reuss-Hill或者Hashin-Shtrikman平均计算混合矿物的弹性模量；
（２）利用DEM模型将孔隙加进系统，并计算干燥岩石骨架的体积模量和剪切模量；
（３）利用Hudson理论和Schoenberg＆Sayers理论，在碳酸盐岩介质各向同性背景中加入裂缝系统，修正计算出的模量；
（４）利用Wood公式将孔隙流体进行混合，计算出混合流体的体积模量；
（ 5 ）计算裂缝型碳酸盐岩的弹性参数和各向异性参数．

1.利用Voigt-Reuss-Hill或者Hashin-Shtrikman平均计算混合矿物的弹性模量

利用岩石基质组分，计算混合矿物的弹性模量

井的矿物成分含量和所需的井曲线如下：

经过Hashin-Shtrikman计算后获得的等效矿物模型结果如下：

//////////////////////////////////////////////////////////////////////////////////////////////
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;

}

２.利用DEM模型将孔隙加进系统，并计算干燥岩石骨架的体积模量和剪切模量；

孔隙空间分成主成分孔隙和次成分孔隙

利用广义xuwhite模型计算岩石的弹性模量：

其中：

采用微分等效介质模型加入孔隙。

\[(1-y)d[K^*(y)]/dy=(K_2-K^*)P^{*2}(y) \]

\[(1-y)d[u^*(y)]/dy=(u_2-u^*)Q^{*2}(y) \]

初始条件是

K(0)=K1和μ（0）=μ1；

其中：

K1，μ1为初始主相材料的体积模量和剪切模量

K2，μ2为逐渐加入包含物材料的体积模量和剪切模量

y=相2的含量

对于流体包含物和空包含物，y等于孔隙度φ，项P和Q是给定的几何因素，P和Q的上标*2指的是几何因素具有等效的背景介质中的包含物材料（具体背景材料参数和系数参考岩石物理手册P110）。

DEM计算后的结果：

​ 紫色为基质正演结果，红色为干骨架正演结果，绿色为实测数据； 
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.利用Ｈｕｄｓｏｎ理论和Ｓｃｈｏｅｎｂｅｒｇ ＆ Ｓａｙｅｒｓ理论，在碳酸盐岩介质各向同性背景中加入裂缝系统，修正计算出的模量。

从前面1,2步建立的等效岩石背景基质中，获得基质的体积模量和剪切模量。在其基础上加入裂缝的影响，修正计算的模量。

目前有两种应用广泛的裂缝模型，Ｈｕｄｓｏｎ 模型和Ｓｃｈｏｅｎｂｅｒｇ线性滑动模型。Ｈｕｄｓｏｎ理论是建立在整体平均的基础上。该模型假设弹性介质的内部分布着薄硬币形状的椭球缝隙。Ｈｕｄｓｏｎ模型利用裂缝密度和缝隙纵横比来描述裂缝系统的结构。裂缝密度e的定义如下：

\[e = \frac{{3\phi }}{{4\pi a}} \]

其中，φ代表孔隙度，α 是裂隙纵横比．则Ｈｕｄｓｏｎ模型的弹性矩阵犆C表示为：

\[C^{eff}_{ij}=C^{0}_{ij}+C^{1}_{ij}+C^{2}_{ij} \]

式中，C0为各向同性背景弹性矩阵，C1,C2为引入裂隙的各向异性一阶校正和二阶校正。

Ｓｃｈｏｅｎｂｅｒｇ模型的弹性矩阵由各向同性参数λ 和μ，以及垂直于裂缝面和平行于裂缝面的弹性差值ΔＮ 和ΔＴ 进行描述。表示裂缝对骨架岩石各向同性特征的影响。在裂缝诱导等效各向异性介质中，裂缝等效介质弹性参数矩阵C****表示为背景各向同性介质弹性参数Cb与裂缝介质弹性参数Cf之差；Ｓｃｈｏｅｎｂｅｒｇ线性滑动模型的弹性矩阵为：

△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公式进行的平均等效模型表达式为：

\[\begin{array}{l} V = \sqrt {\frac{{{K_R}}}{\rho }} \\ \frac{1}{{{K_R}}} = \sum\limits_{i = 1}^N {\frac{{{f_i}}}{{{K_i}}}} \\ \rho = \sum\limits_{i = 1}^N {{f_i}} {\rho _i} \end{array}\ \]

其中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 ）计算裂缝型碳酸盐岩的弹性参数和各向异性参数．

由于研究区域只含有气，不用进行流体替换。

各向异性参数的表达式为：

\[\begin{array}{l} \varepsilon = \left( {C[1][1] - C[3][3]} \right)/\left( {2*C[3][3]} \right);\\ \delta = {\rm{((C[1][3] + C[5][5])*(C[1][3] + C[5][5]) - (C[3][3] - C[5][5])*(C[3][3] - C[5][5])) / (2 * C[3][3] * (C[3][3] - C[5][5]));}}\\ \gamma {\rm{ = (C[6][6] - C[4][4]) / (2 * C[4][4]);}} \end{array} \]

对应的计算代码：

		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]);

计算所得的各向异性参数与实测的吻合度较高。

参考及其引用文献：

张广智，陈怀震，王琪等．基于碳酸盐岩裂缝岩石物理模型的横波速度和各向异性参数预测．地球物理学报，２０１３，５６（５）：１７０７
１７１５，ｄｏｉ：１０．６０３８／ｃｊｇ２０１３０５２８．

岩石物理手册.Gray Mavko and so on.

标签：CHTI,流程,float,Por,HS2,HS1,碳酸盐岩,裂缝
From： https://www.cnblogs.com/GeophysicsWorker/p/17046950.html

相关文章

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99