cementedsand.m

来自「Stanford的SRB实验室Quantitative Seismic Inte」· M 代码 · 共 141 行

function [Phic,Rhodc,Vpdc,Vsdc,Rhosc,Vpsc,Vssc,Phis,Rhods,Vpds,Vsds,Rhoss,Vpss,Vsss]=cementedsand(Clay,Feldspar,Calcite,phic,fcement,kf,rhof,shearfact)
%        [Phic,Rhodc,Vpdc,Vsdc,Rhosc,Vpsc,Vssc,Phis,Rhods,Vpds,Vsds,Rhoss,Vpss,Vsss]=cementedsand(Clay,Feldspar,Calcite,phic,fcement,kf,rhof,shearfact)
%
% Calculates the dry- and saturated rock M-modulus, shear modulus and
% Poisson's ratio.  Model uses contact cement theory to set the
% right end member starting from critical porosity, and the  modified upper Hashin-Shtrikman 
% model to continue cement model to lower porosities.  Finally, a line of
% constant cement is generated with modified lower bound, beginning from a
% point on the cement line.
% The zero-porosity end point is pure mineral mix.
%
% INPUTS
%   Clay			Volume clay content in solid phase (fraction)
%   Feldspar		Volume feldspar content in solid phase (fraction)
%   Calcite		    Volume calcite content in solid phase (fraction)
%   PhiC			Critical porosity ~0.4
%   fcement         percent fraction
%   Kf			    Pore fluid bulk modulus
%   RHOf			Pore fluid density
%   shearfact       empirical correction to contact cement shear modulus
%
% OUTPUTS
%   PhiC		    Porosity (fraction) -- the same as in inputs
%   Rhodc           Dry rock bulk density, cemented sand
%   Vpdc            Dry rock Vp, cemented sand
%   Vsdc            Dry rock Vs, cemented sand
%   Rhosc           Saturated rock bulk density, cemented sand
%   Vpsc            Saturated rock Vp, cemented sand
%   Vssc            Saturated rock Vs, cemented sand
%   Rhods           Dry rock bulk density, soft shale
%   Vpds            Dry rock Vp, soft shale
%   Vsds            Dry rock Vs, soft shale
%   Rhoss           Saturated rock bulk density, soft shale
%   Vpss            Saturated rock Vp, soft shale
%   Vsss            Saturated rock Vs, soft shale
%Quartz content is 1-Clay-Feldspar-Calcite
%Written by Gary Mavko
% Mineral component elastic moduli
KQuartz=36.6e9; GQuartz=45e9; KClay=25e9; GClay=9e9; KFeldspar=75.6e9; GFeldspar=25.6e9;
KCalcite=76.8e9; GCalcite=32.0e9;
% Mineral component densities
RhoQuartz=2650; RhoClay=2580; RhoFeldspar=2630; RhoCalcite=2710;
if nargin<1,
    prompt = {'critical porosity',...
             'cement fraction', ...
             'Clay Fraction', ...
             'Calcite Fraction', ...
             'Feldspar Fraction', ...
             'Shear Reduction Factor', ...
             'Fluid Bulk Modulus',...
             'Fluid Density'};    
    defans = {'.4', ...
              '.03', ...
             '.05',...
             '0',...
             '0', ...
             '.5', ...
             '2.5e9',...
             '1020'};
    answer    = inputdlg(prompt,'Soft Sediment Model',1,defans);
    phic      = str2num(answer{1});
    fcement   = str2num(answer{2});
    Clay      = str2num(answer{3});
    Calcite   = str2num(answer{4});
    Feldspar  = str2num(answer{5});
    shearfact = str2num(answer{6});
    kf        = str2num(answer{7});
    rhof      = str2num(answer{8});
end;
% Mineral phase moduli from Voigt-Reuss-Hill average
Quartz=1-Clay-Feldspar-Calcite;
kmin = 0.5.*(Quartz.*KQuartz+Clay.*KClay+Feldspar.*KFeldspar+Calcite.*KCalcite+1./(Quartz./KQuartz+Clay./KClay+Feldspar./KFeldspar+Calcite./KCalcite));
gmin = 0.5.*(Quartz.*GQuartz+Clay.*GClay+Feldspar.*GFeldspar+Calcite.*GCalcite+1./(Quartz./GQuartz+Clay./GClay+Feldspar./GFeldspar+Calcite./GCalcite));
mmin = kmin + (4/3)*gmin;
% Mineral phase density
rhomin=Quartz.*RhoQuartz+Clay.*RhoClay+Feldspar.*RhoFeldspar+Calcite.*RhoCalcite;
% Compute Saturated Hashin-Shtrikman  lower bound at PhiC.  Lower is the suspension line
%[kuc,klc,guc,glc]=hsbound(kmin,gmin,kf,0,phic);
% create modified diagenetic trend, which is a combination of the Dvorkin cement model
% and modified upper Hashin-Shtrikman bound.
[kcemsat, gcemsat, rhocemsat, phicem] = diagen(phic, kmin, gmin, rhomin, kf, rhof, shearfact);
mcemsat  = kcemsat + (4/3).*gcemsat;  
Vpsc = sqrt(mcemsat./rhocemsat);      
Vssc = sqrt(gcemsat./rhocemsat);     
[kcemdry, gcemdry, rhocemdry, phicem] = diagen(phic, kmin, gmin, rhomin, 0, 0, shearfact);
mcemdry  = kcemdry + (4/3).*gcemdry; 
Vpdc     = sqrt(mcemdry./rhocemdry); 
Vsdc     = sqrt(gcemdry./rhocemdry); 
% modified lower HS bounds for sorting lines
    phis      = phic - fcement;
    indsor    = (phicem<=phis);
    phisor    = phicem(indsor);
    rhosordry = rhocemdry(indsor);
    rhosorsat = rhocemsat(indsor);
    kssat     = min(kcemsat(indsor));
    ksdry     = min(kcemdry(indsor));
    gssat     = min(gcemsat(indsor));
    gsdry     = min(gcemdry(indsor));
    mssat     = kssat + (4/3)*gssat;
    msdry     = ksdry + (4/3)*gsdry;
    % Lower HS
    khssat=1./((phisor./phis)./(kssat+4.*gssat./3)+((phis-phisor)./phis)./(kmin+4.*gssat./3))-4.*gssat./3;
    khsdry=1./((phisor./phis)./(ksdry+4.*gsdry./3)+((phis-phisor)./phis)./(kmin+4.*gsdry./3))-4.*gsdry./3;
    ZZ1=(gssat./6).*(9.*kssat+8.*gssat)./(kssat+2.*gssat);
    ghssat=1./((phisor./phis)./(gssat+ZZ1)+((phis-phisor)./phis)./(gmin+ZZ1))-ZZ1; 
    ZZ1=(gsdry./6).*(9.*gsdry+8.*gsdry)./(ksdry+2.*gsdry);
    ghsdry=1./((phisor./phis)./(gsdry+ZZ1)+((phis-phisor)./phis)./(gmin+ZZ1))-ZZ1; 
    Vpss  = sqrt((khssat+(4/3)*ghssat)./rhosorsat);
    Vpds  = sqrt((khsdry+(4/3)*ghsdry)./rhosordry);
    Vsss  = sqrt(ghssat./rhosorsat);
    Vsds  = sqrt(ghsdry./rhosordry);
    j = phicem >= (phic - 2*fcement);
    phicem  = phicem(j);
    kcemsat = kcemsat(j);
    gcemsat = gcemsat(j);
    rhocemsat = rhocemsat(j);
    mcemsat  = kcemsat + (4/3).*gcemsat;  
    kcemdry = kcemdry(j);
    gcemdry = gcemdry(j);
    rhocemdry = rhocemdry(j);
    mcemdry  = kcemdry + (4/3).*gcemdry;  
    Vpsc = sqrt(mcemsat./rhocemsat);      
    Vssc = sqrt(gcemsat./rhocemsat);     
    mcemdry  = kcemdry + (4/3).*gcemdry; 
    Vpdc     = sqrt(mcemdry./rhocemdry); 
    Vsdc     = sqrt(gcemdry./rhocemdry); 
    Phic   = phicem;
    Phis   = phisor;
    Rhodc = rhocemdry;
    Rhosc = rhocemsat;
    Rhods = rhosordry;
    Rhoss = rhosorsat;
% -----------------------------------------------------
function [ku,kl,gu,gl] = hsbound(k1,mu1,k2,mu2,por)
	ku=k2+(1.-por)*(k1-k2)*(k2+4.*mu1/3.)./(k2+4.*(mu1/3.)+por*(k1-k2));
	kl=k2+(1.-por)*(k1-k2)*(k2+4.*mu2/3.)./(k2+4.*(mu2/3.)+por*(k1-k2));
	fgu=mu1*(9.*k1+8.*mu1)/(6.*(k1+2.*mu1));
	fgl=mu2*(9.*k2+8.*mu2)/(6.*(k2+2.*mu2));
	gu=mu2+(mu1-mu2)*(1.-por)*(mu2+fgu)./(mu2+fgu+por*(mu1-mu2));
	gl=mu2+(mu1-mu2)*(1.-por)*(mu2+fgl)./(mu2+fgl+por*(mu1-mu2));