function [Vp1,Vp2,Vs,Qp1,Qp2,Qs,w,Sw]=BISQ(K,Ks,Kw,Kg,de_w,de_s,de_a,de_g,np,k,R,v,f)
% 该函数利用BISQ模型计算快纵波、慢纵波和横波的速度随孔隙度曲线，快纵波和慢纵波品质因子曲线。
% 输出：
% Vp1-快纵波随孔隙度曲线
% Vp2-慢纵波随孔隙度曲线
% Vs-横波速度随孔隙度曲线
% Qp1-快纵波品质因子倒数曲线
% Qp2-慢纵波品质因子倒数曲线
% Qs-横波逆品子因子曲线
% 编写人：朱帅润 成都理工大学
format long
K=16*10^9;
Ks=38*10^9; % Ks-岩石颗粒的体积模量
Kw=2.19*10^9; % Kw-空隙流体水的体积模量
Kg=1.42*10^5; % Kg-空隙流体气体的体积模量
de_w=1000; % density_f-空隙中流体水密度
de_g=140; % density_f-空隙中空气密度
de_s=2650; % density_s-固体骨架的密度
de_a=420; % density_a-附加质量的密度
v=0.15; % 泊松比
np=0.001; % rp-流体粘度
k=1.25*10^-15; % k-渗透率
R=0.001; % R-特征喷射流动长度单位m
f=1552; % 频率取值Hz
Opro=0.01:0.01:0.4; %孔隙度取值
Sw=0.01:0.18:0.91; %饱和度赋值
in=length(Opro);
hn=length(Sw);
K_s=zeros(hn,in); %饱和体积模量
K_d=zeros(hn,in); %干燥体积模量
lamu=zeros(hn,in);
S_q=zeros(hn,in);
F=zeros(hn,in);
Seta=zeros(hn,in);
A=zeros(hn,in);
C=zeros(hn,in);
B=zeros(hn,in);
P_x=zeros(hn,in);
M=zeros(hn,in);
G=zeros(hn,in);
X1=zeros(hn,in);
X2=zeros(hn,in);
Y=zeros(hn,in);
Vp1=zeros(hn,in);
Vp2=zeros(hn,in);
Vs=zeros(hn,in);
Qp1=zeros(hn,in);
Qp2=zeros(hn,in);
Qs=zeros(hn,in);
for hh=1:hn
w=2*pi*f; %角频率
de_f=Sw(hh)*de_w+(1-Sw(hh))*de_g;
K_f=Sw(hh)*Kw+(1-Sw(hh))*Kg;
Rs=R*sqrt(Sw(hh));
for ii=1:in
a=1-K/Ks;
M(hh,ii)=1/(Opro(ii)*0.25/K_f+(a-Opro(ii))/Ks);
G(hh,ii)=0.5*M(hh,ii)*(1-2*v)/(1-v);
%G(hh,ii)=70.369*Opro(ii)*Opro(ii)-55.312*Opro(ii)+12.981;
%M(hh,ii)=2.0*G(hh,ii)*(1-v)/(1-2*v);
% K_s(hh,ii)=59.044*Opro(ii)*Opro(ii)-65.499*Opro(ii)+25.718;
% TEMPone=1/(K_s(hh,ii)/(Ks-K_s(hh,ii))-K_f/(Opro(ii)*(Ks-K_s(hh,ii))))+1;
% K_d(hh,ii)=1/TEMPone*Ks;
temp1=1/K_f+(a-Opro(ii))/(Opro(ii)*Ks);
F(hh,ii)=1/temp1;
temp2=(de_a/de_f+Opro(ii))/Opro(ii)+1i*(Opro(ii)*np/(w*de_f*k));
Seta(hh,ii)=1/temp2;
temp3=de_f*w*w/(F(hh,ii)*Seta(hh,ii));
lamu(hh,ii)=sqrt(temp3);
temp=1-2.0*besselj(1,R*lamu(hh,ii))/(lamu(hh,ii)*R*besselj(0,R*lamu(hh,ii)));
S_q(hh,ii)=temp*Sw(hh);
P_x(hh,ii)=(1-Opro(ii))*de_s+(1-Seta(hh,ii))*Opro(ii)*de_f;
C(hh,ii)=((1-Opro(ii))*de_s/(Opro(ii)*de_f))+(1+((1-Opro(ii))*de_s/(Opro(ii)*de_f)))* ...
(de_a/(Opro(ii)*de_f)+1i*((np*Opro(ii))/(k*de_f*w)));
temp4=F(hh,ii)*S_q(hh,ii)*(2*a-Opro(ii)-Opro(ii)*((1-Opro(ii))*de_s/(Opro(ii)*de_f)))- ...
(M(hh,ii)+F(hh,ii)*S_q(hh,ii)*a*a/Opro(ii))*(1+de_a/(Opro(ii)*de_f)+1i*((np*Opro(ii))/(k*de_f*w)));
B(hh,ii)=temp4/(Opro(ii)*de_f);
A(hh,ii)=Opro(ii)*F(hh,ii)*S_q(hh,ii)*M(hh,ii)/(Opro(ii)*de_f*Opro(ii)*de_f);
X1(hh,ii)=sqrt(-B(hh,ii)/(2*A(hh,ii))+sqrt((B(hh,ii)/(2*A(hh,ii)))^2-C(hh,ii)/A(hh,ii)));
X2(hh,ii)=sqrt(-B(hh,ii)/(2*A(hh,ii))-sqrt((B(hh,ii)/(2*A(hh,ii)))^2-C(hh,ii)/A(hh,ii)));
Y(hh,ii)=sqrt(P_x(hh,ii));
Vp1(hh,ii)=1/real(X1(hh,ii));
Vp2(hh,ii)=1/real(X2(hh,ii));
Vs(hh,ii)=sqrt(G(hh,ii))/real(Y(hh,ii));
Qp1(hh,ii)=2*imag(X1(hh,ii))/real(X1(hh,ii));
Qp2(hh,ii)=2*imag(X2(hh,ii))/real(X2(hh,ii));
Qs(hh,ii)=2*imag(Y(hh,ii))/real(Y(hh,ii));
end
end
figure(1)
plot(Opro*100,Vp1(2,:),Opro*100,Vp1(3,:),Opro*100,Vp1(4,:),Opro*100,Vp1(5,:),Opro*100,Vp1(6,:),'Linewidth',1.5);
grid on;
xlabel('孔隙度%');
ylabel('快Vp m/s');
legend('S=19%','S=37%','S=55%','S=73%','S=91%');
figure(2)
plot(Opro*100,Vp2(1,:),Opro*100,Vp2(2,:),Opro*100,Vp2(3,:),Opro*100,Vp2(4,:),Opro*100,Vp2(5,:),Opro*100,Vp2(6,:),'Linewidth',1.5);
grid on;
xlabel('孔隙度%');
ylabel('慢Vp m/s');
legend('S=1%','S=19%','S=37%','S=55%','S=73%','S=91%');
figure(3)
plot(Opro*100,Vs(1,:),Opro*100,Vs(2,:),Opro*100,Vs(3,:),Opro*100,Vs(4,:),Opro*100,Vs(5,:),Opro*100,Vs(6,:),'Linewidth',1.5);
grid on;
xlabel('孔隙度%');
ylabel('Vs m/s');
legend('S=1%','S=19%','S=37%','S=55%','S=73%','S=91%');
figure(4)
plot(Sw,Qp1(:,10),Sw,Qp1(:,15),Sw,Qp1(:,20),Sw,Qp1(:,25),Sw,Qp1(:,30),Sw,Qp1(:,38),'Linewidth',1.5);
grid on;
xlabel('饱和度');
ylabel('1/Q');
legend('孔隙度=10%','孔隙度=15%','孔隙度=20%','孔隙度=25%','孔隙度=30%','孔隙度=38%');
figure(5)
plot(Opro*100,Qp1(1,:),Opro*100,Qp1(2,:),Opro*100,Qp1(3,:),Opro*100,Qp1(4,:),Opro*100,Qp1(5,:),Opro*100,Qp1(6,:),'Linewidth',1.5);
grid on;
xlabel('孔隙度%');
ylabel('1/Q');
legend('S=1%','S=19%','S=37%','S=55%','S=73%','S=91%');
figure(6)
plot(Sw,Qp2(:,10),Sw,Qp2(:,15),Sw,Qp2(:,20),Sw,Qp2(:,25),Sw,Qp2(:,30),Sw,Qp2(:,38),'Linewidth',1.5);
grid on;
xlabel('饱和度');
ylabel('1/Q');
legend('孔隙度=10%','孔隙度=15%','孔隙度=20%','孔隙度=25%','孔隙度=30%','孔隙度=38%');
figure(7)
plot(Sw,Qs(:,10),Sw,Qs(:,15),Sw,Qs(:,20),Sw,Qs(:,25),Sw,Qs(:,30),Sw,Qs(:,38),'Linewidth',1.5);
grid on;
xlabel('饱和度');
ylabel('1/Q');
legend('孔隙度=10%','孔隙度=15%','孔隙度=20%','孔隙度=25%','孔隙度=30%','孔隙度=38%');
评论0