I'm writing some functions in numpy for rock physics modelling and have noticed that one of my functions gives erroneous results. The function is my implimentation of Hertz-Mindlin sphere modelling:
Summary of the Hertz-Mindlin model
Here is my function currently:
# Hertz-Mindlin sphere pack model:
import numpy as np
def hertzmindlin(K0, G0, PHIC, P, f=1.0):
'''
Hertz-Mindlin sphere-pack model, adapted from:
'Dvorkin, J. and Nur, A., 1996. Elasticity of high-porosity sandstones:
Theory for two North Sea data sets. Geophysics, 61(5), pp.1363-1370."
Arguments:
K0 = Bulk modulus of mineral in GPa
G0 = Shear modulus of mineral in GPa
PHIC = Critical porosity for mineral-fluid mixture. Calculate using Dvorkin-Nuir (1995) or use literature
P = Confining pressure in GPa
f = Shear modulus correction factor. Default = 1
Results:
V0 = Theoretical poissons ratio of mineral
n = Coordination number of sphere-pack, calculated from Murphy's (1982) empirical relation
K_HM = Hertz-Mindlin effective dry Bulk modulus at pressure, P, in GPa
G_HM = Hertz-Mindlin effective dry Shear modulus at pressure, P, in GPa
'''
V0 = (3*K0-2*G0)/(6*K0+2*G0) # Calculated theoretical poissons ratio of bulk rock
n = 20-(34*PHIC)+(14*(PHIC**2)) # Coordination number at critical porosity (Murphy 1982)
K_HM = (P*(n**2*(1-PHIC)**2*G0**2) / (18*np.pi**2*(1-V0)**2))**(1/3)
G_HM = ((2+3*f-V0*(1+3*f))/(5*(2-V0))) * ((P*(3*n**2*(1-PHIC)**2*G0**2)/(2*np.pi**2*(1-V0)**2)))**(1/3)
return K_HM, G_HM
The problem is that when I run this function for inputs of:
K, G, = 36, 45
PHIC = 0.4
P = 0.001
I get a result of K_HM = 1.0, G_HM = 0.49009009009009
The hand calculated and excel calculated values show this is incorrect, I should be outputting K_HM = 0.763265313, G_HM = 1.081083984
I am fairly certain something is going wrong in the function based on the fact that for the inputs K, G, the output G should be larger than K (it is currently smaller)
Any help would be appreciated! I can do this in excel, but ideally want everything running in python.
In Python2, division of integers (using /) returns an integer. For example, 1/3 = 0.
In Python3, division of integers (using /) may return a float.
It appears you are using Python2. To get floating-point division (in both Python2 and Python3), ensure each division operation involves at least one float: for example, change 1/3 to 1.0/3 or 1/3.0 or (acceptable but perhaps less readable, 1/3.):
import numpy as np
def hertzmindlin(K0, G0, PHIC, P, f=1.0):
K0, G0 = map(float, (K0, G0))
V0 = (3*K0-2*G0)/(6*K0+2*G0) # Calculated theoretical poissons ratio of bulk rock
n = 20-(34*PHIC)+(14*(PHIC**2)) # Coordination number at critical porosity (Murphy 1982)
K_HM = (P*(n**2*(1-PHIC)**2*G0**2) / (18*np.pi**2*(1-V0)**2))**(1/3.0)
G_HM = ((2+3*f-V0*(1+3*f))/(5*(2-V0))) * ((P*(3*n**2*(1-PHIC)**2*G0**2)/(2*np.pi**2*(1-V0)**2)))**(1/3.0)
return K_HM, G_HM
K, G, = 36, 45
PHIC = 0.4
P = 0.001
print(hertzmindlin(K, G, PHIC, P))
Alternatively, in later versions of Python2 (such as Python2.7) you could place
from __future__ import division
at the top of your script (before all other import statements) to activate Python3-style floating-point division.