2017-12-05 6 views
2

나는 바위 물리 모델링을 위해 numpy에서 몇 가지 함수를 작성하고 있으며 내 함수 중 하나가 잘못된 결과를 제공한다는 것을 알고 있습니다. 여기 잘못된 결과를주는 numpy 함수 - 손으로 확인하고 Excel을

Summary of the Hertz-Mindlin model

현재 내 함수 : 함수는 헤르츠-Mindlin 구형 모델의 내도 구현입니다

# 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 

문제는 그 난의 입력에 대해이 기능을 실행하면

K, G = 36, 45

PHIC = 0.4

,515,

P = 0.001

난 K_HM = 1.0 G_HM 결과를 얻을 = 0.49009009009009 손이 계산되고 계산 된 값이 잘못 표시 엑셀

, I는 K_HM = 0.763265313, G_HM = 1.081083984

출력되어야

나는 입력 K, G에 대해 출력 G가 K보다 커야한다는 사실에 근거하여 함수에서 잘못된 것이 확실하다. (현재는 더 작다)

어떤 도움을 주실 수 있겠습니까? Excel에서이 작업을 수행 할 수 있지만, 모든 것이 Python으로 실행되기를 원합니다.

답변

2

파이썬 2에서 정수의 나눗셈 (/ 사용)은 정수를 반환합니다. 예 : 1/3 = 0. Python3에서 정수의 나눗셈 (/ 사용)은 실수를 반환 할 수 있습니다.

Python2를 사용하고있는 것으로 보입니다. 예 1/3-1.0/3 또는 1/3.0 또는 (허용하지만 아마도 덜 판독 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)) 
(Python2 및 Python3 모두) 부동 소수점 나눗셈을 얻기 위해, 각각의 분할 동작이 적어도 하나 개의 유동을 포함 확보

또는 (예 : Python2.7 등) Python2의 이후 버전에서 당신은에 (다른 모든 import 문 전) 스크립트의 상단에

from __future__ import division 

을 배치 할 수 있습니다.

+0

굉장! 그것은 숫자 형식과 관련된 것이라고 생각합니다. 불행히도 내 교육 기관에서 내 파이썬 버전이나 환경을 변경할 수 없지만 나중에 이것을 고려할 것입니다! –