Source code for bruges.rockphysics.rpt

Rock physics templates
Uses rock physics models defined in
For Voigt-Reuss-Hill averages I define my own function here but I should
really change it and use the bounds already defined in bruges.
@author: Alessandro Amato del Monte
import numpy as np
import matplotlib.pyplot as plt
from bruges.rockphysics import rpm
from bruges.rockphysics import fluidsub

[docs]def vrh(volumes,k,mu): ''' Calculates Voigt-Reuss-Hill bounds. Args: volumes: array with volumetric fractions k: array with bulk modulus mu: array with shear modulus Returns: k_u, k_l: upper (Voigt) and lower (Reuss) average for k mu_u, mu_l: upper (Voigt) and lower (Reuss) average for mu k0, mu0: Hill average of k and mu ''' f=np.array(volumes).T k=np.resize(np.array(k),np.shape(f)) mu=np.resize(np.array(mu),np.shape(f)) ax=0 if f.ndim==1 else 1 k_u = np.sum(f*k,axis=ax) k_l = 1./np.sum(f/k,axis=ax) mu_u = np.sum(f*mu,axis=ax) mu_l = 1./np.sum(f/mu,axis=ax) k0 = (k_u+k_l)/2. mu0 = (mu_u+mu_l)/2. return k_u, k_l, mu_u, mu_l, k0, mu0
[docs]def rpt(model='soft',vsh=0.0,fluid='gas',phic=0.4,Cn=8,P=10,f=1,cement='quartz'): if cement=='quartz': Kc, Gc = 37, 45 elif cement=='calcite': Kc, Gc = 76.8, 32 elif cement=='clay': Kc, Gc = 21, 7 phi=np.linspace(0.1,phic-.1,6) sw=np.linspace(0,1,5) (Khc, Dhc) = (Kg, Dg) if fluid == 'gas' else (Ko,Do) K0,G0 = vrh([vsh, 1-vsh],[Ksh,Kqz],[Gsh,Gqz])[4:] D0 = vsh*Dsh+(1-vsh)*Dqz if model=='soft': Kdry, Gdry = rpm.softsand(K0,G0,phi,phic,Cn,P,f) elif model=='stiff': Kdry, Gdry = rpm.stiffsand(K0,G0,phi,phic,Cn,P,f) elif model=='cem': Kdry, Gdry = rpm.contactcement(K0,G0,phi,phic,Cn,Kc,Gc,scheme=2) elif model=='crit': Kdry, Gdry = rpm.critpor(K0,G0,phi,phic) xx=np.empty((phi.size,sw.size)) yy=np.empty((phi.size,sw.size)) for i,val in enumerate(sw): Kf = vrh([val,1-val],[Kb,Khc],[999,999])[1] Df = val*Db+(1-val)*Dhc vp,vs,rho,_= fluidsub.vels(Kdry,Gdry,K0,D0,Kf,Df,phi) xx[:,i]=vp*rho yy[:,i]=vp/vs plt.figure(figsize=(10,6)) plt.plot(xx, yy, '-ok', alpha=0.3) plt.plot(xx.T, yy.T, '-ok', alpha=0.3) for i,val in enumerate(phi): plt.text(xx[i,-1],yy[i,-1]+.01,'$\phi={:.02f}$'.format(val), backgroundcolor='0.9') plt.text(xx[-1,0]-100,yy[-1,0],'$S_w={:.02f}$'.format(sw[0]),ha='right', backgroundcolor='0.9') plt.text(xx[-1,-1]-100,yy[-1,-1],'$S_w={:.02f}$'.format(sw[-1]),ha='right', backgroundcolor='0.9') plt.xlabel('Ip'), plt.ylabel('Vp/Vs') plt.xlim(xx.min()-xx.min()*.1,xx.max()+xx.max()*.1) plt.ylim(yy.min()-yy.min()*.1,yy.max()+yy.max()*.1) plt.title('RPT (N:G={0}, fluid={1})'.format(1-vsh, fluid))