cement.m

Stanford的SRB实验室Quantitative Seismic Interpretation的免费MATLAB程序

function [ksat, gframe] = cement(varargin)
%CEMENT Cemented sand model (different inputs)
%Calculates the water-saturated elastic moduli versus porosity lines
%using Dvorkin's contact cement model
%Inputs are by input dialog box, if called without input arguments.
%     phi:  porosity values where moduli are to be computed
%     phic: critical porosity
%     c:    coordination number
%     gs:   grain shear modulus
%     ks:   grain bulk modulus
%     gc:   cement shear modulus
%     kc:   cement bulk modulus
%     kf:   fluid bulk modulus
%     scheme: (1 or 2) Cementation scheme 1: cement at contact, 2: on surface
%Plots effective moduli vs. porosity when called with no output arguments.


%Written by Jack Dvorkin
%I/O modifications, Tapan Mukerji, 6/1999.
%ref. Dvorkin & Nur, 1996, Geophysics, 61, 1363-1370
%     Rock Physics Handbook, section 5.2
%modified by Gary Mavko 3/01 to input porosities for calculation and change I/O to K, G

prompt={'Critical Porosity','Coord.#','G-grain','K-grain','G-cement','K-cement','K-fluid','Cem. Scheme (1 or 2)' };
defans={'.38','8.5','45e9','36e9','45e9','36e9','2.2e9','2'};

Phi0 = varargin{1};
if nargin<2
    getpar=inputdlg(prompt,'Contact Cement Model',1,defans);
    for k=1:length(getpar), param(k)=str2num(getpar{k}); end;
    phic=param(1); c=param(2); g=param(3); k=param(4); gc=param(5); kc=param(6);
    kf=param(7); schopt=param(8);
else
    phic=varargin{2}; c=varargin{3}; g=varargin{4}; k=varargin{5};
    gc=varargin{6}; kc=varargin{7}; kf=varargin{8}; schopt=varargin{9};
end;

nu = 0.5*(3*k-2*g)./(3*k+g);
nuc = 0.5*(3*kc-2*gc)./(3*kc+gc);

% Fraction of cement in the rock
       fc = phic-Phi0;
% Fraction of grain in the solid
       fgs = (1-phic)./(1-Phi0);
% Fraction of cement in the solid
       fcs = (phic-Phi0)./(1-Phi0);
% Bulk modulus of the solid;   Voigt - Reuss - Hill average
       ks = (fgs.*k+fcs.*kc+1./(fgs./k+fcs./kc))./2;
% Shear modulus of the solid
       gs = (fgs.*g+fcs.*gc+1./(fgs./g+fcs./gc))./2;
% Kframe and Gframe at Phi=Phi0, Cementation Scheme 1
       if schopt==1
           a = 2.*(((phic-Phi0)./(3.*c.*(1.-phic))).^0.25);
       end;
% Cementation Scheme 2
       if schopt==2
           a = sqrt((2.*(phic-Phi0))./(3.*(1.-phic)));
       end;
% Capital Lambdas
       alam = (2./pi).*(gc./g).*(1.-nu).*(1.-nuc)./(1.-2.*nuc);
       alamtau = (1./3.14).*(gc./g);
% Effective bulk modulus
       r1 = kc+4.*gc./3;
       r2 = c.*(1.-phic)./6;
       r3 = -0.024153.*(alam.^(-1.3646)).*(a.^2)+0.20405.*(alam.^(-0.89008)).*a+0.00024649.*(alam.^(-1.9864));
       kframe = r1.*r2.*r3;
% Effective shear modulus
       r1 = gc;
       r2 = 3.*c.*(1.-phic)./20;
       a1t = -0.01.*(2.2606.*nu.*nu+2.0696.*nu+2.2952);
       a2t = 0.079011.*nu.*nu+0.17539.*nu-1.3418;
       b1t = 0.05728.*nu.*nu+0.09367.*nu+0.20162;
       b2t = 0.027425.*nu.*nu+0.052859.*nu-0.87653;
       c1t = 0.0001.*(9.6544.*nu.*nu+4.9445.*nu+3.1008);
       c2t = 0.018667.*nu.*nu+0.4011.*nu-1.8186;
       r3 = (a1t.*(alamtau.^a2t)).*a.*a+(b1t.*(alamtau.^b2t)).*a+c1t.*(alamtau.^c2t);
       gframe = 0.6.*kframe + r1.*r2.*r3;

	  kframe(end) = 0;
	  gframe(end) = 0;
% Gassmann
       ksat = ks.*(Phi0.*kframe-(1+Phi0).*kf.*kframe./ks+kf)./((1-Phi0).*kf+Phi0.*ks-kf.*kframe./ks); 


if nargout==0
     subplot(1,2,1)
     plot(Phi0,ksat,'r-')
     set(gca,'fontname','bookman','fontsize',9)
     xlabel('Porosity','fontname','bookman','fontsize',11)
     ylabel('Bulk Modulus','fontname','bookman','fontsize',11)
     subplot(1,2,2)
     plot(Phi0,gframe,'r-')
     set(gca,'fontname','bookman','fontsize',9)
     xlabel('Porosity','fontname','bookman','fontsize',11)
     ylabel('Shear Modulus','fontname','bookman','fontsize',11)
end;