#include <stdio.h>
#include <math.h>
int main ()
{
double Km,um,a,b,c,d,Pi,Qi,Kkt,ukt,fs,Sw,Kf,uf,densf,densma,dens,Ksat,usat,Vp,Vs;
double pai=3.14159;
int i,j;
scanf("%lf%lf%lf%lf",&fs,&c,&Sw,&i);//需知砂岩的体积比例，孔隙度,水的饱和度，求干岩石令Xi=Ki=ui=0
for(j=1;j<=100;j++)
{
Km=(1.0/2)*(37*fs+20.9*(1-fs)+1.0/(fs/37+(1-fs)/20.9)); //利用V-R-H求取骨架体积模量，剪切模量
um=(1.0/2)*(44*fs+6.85*(1-fs)+1.0/(fs/44+(1-fs)/6.85));
d=(um/6)*(9*Km+8*um)/(Km+2*um);
a=um*(3*Km+um)/(3*Km+4*um);
b=um*(3*Km+um)/(3*Km+7*um);
Kf=1.0/(Sw/2.2+(1-Sw)/2.5);//已知水的饱和度Sw，利用WOOD公式求取混合流体的弹性模量和密度，基质的密度
uf=0;
densf=Sw+(1-Sw)*1.03; //混合液体密度
densma=2.65*fs+(1-fs)*2.05; //岩石基质密度
dens=(1-c)*densma+c*densf; //饱含孔隙流体岩石的等效密度
{
if (i<2) //spheres
{
Pi=(Km+(4.0/3)*um)/(Kf+(4.0/3)*um);
Qi=(um+d)/(uf+d);
}
else
{
if(i<3) //needles
{
Pi=(Km+um+(1.0/3)*uf)/(Kf+um+(1.0/3)*uf);
Qi=(1.0/5)*(4.0*um/(um+uf)+2*(um+b)/(uf+b)+(Kf+(4.0/3)*um)/(Kf+um+(1.0/3)*uf));
}
else
{
if(i<4) //disks
{
Pi=(Km+(4.0/3)*uf)/(Kf+(4.0/3)*uf);
Qi=(um+d)/(uf+d);
}
else // penny cracks
Pi=(Km+(4.0/3)*uf)/(Kf+(4.0/3)*uf+pai*a*0.01);
Qi=(1.0/5)*(1+8*um/(4*uf+pai*0.01*(um+2*a))+2*(Kf+(2.0/3)*(uf+um))/(Kf+(4.0/3)*uf+pai*0.01*a));
}
}
}
Kkt=Km*pow(1-c,Pi); //KT公式求取干岩石弹性模量
ukt=um*pow(1-c,Qi);
Ksat=Km*(Km/(Km-Kkt)+Kf/((Km-Kf)*c))/(1+Km/(Km-Kkt)+Kf/((Km-Kf)*c));//利用Gassmann求解饱含流体岩石的等效弹性模量
usat=ukt;
Vp=sqrt((Ksat+(4.0/3)*usat)/dens);
Vs=sqrt(usat/dens);
printf("Vp=%f,Vs=%f\n",Vp,Vs);
fs=fs+0.01;
}
return 0;
}
求大神