new

import numpy as np
import matplotlib.pyplot as plt
import math as m
#Definiendo par'ametros de entrada
N=4.016
Lambda=0.348
M=0.941
nu=0.3
Gamma=3.720
Kappa=0.052
H=0.678
p0=300  #kPa#
q0=0    #kPa#
v=1.946
n=10
pf=m.exp((Gamma-v)/Lambda)
qf=M*pf
pfe=(m.exp((Gamma-v)/(Lambda-Kappa)))/(p0**(Kappa/(Lambda-Kappa)))
qfe=M*pfe
pj=m.e*pfe
pe=p0 #por condiciones no drenadas
qe=M*pe*(1-m.log(pe/pfe))
Delta_p=pf-pe
dp=Delta_p/n  #Incremento en p
#-----------
#Pared 1
K=(v*p0)/Kappa
E=3*K*(1-2*nu)
G=E/(2*(1+nu)) #;print(G)
p_i=p0
q_i=qe
j=2
p=[p0]
q=[q0]
d_evole=[0]
for x in range (j,n+2):
    p_f=p_i+(j-1)*dp
    q_f=((M*p_f)/(Lambda-Kappa))*((Lambda-Kappa)-Lambda*m.log(p_f)+(Gamma-v))
    dq=q_f-q_i
    dEvol_e=dp/K
    dEvol_p=((Lambda-Kappa)/(v*M*p_f))*((M-(q_f/p_f))*dp + dq)
    dEvol=dEvol_e+dEvol_p
    dEs_e=dq/(3*G)#; print(dEs_e)
    dEs_p=((Lambda-Kappa)/(v*M*p_f))*(dp+(dq/(M-(q_f/p_f))))
    dEs=dEs_e+dEs_p
    #dEa=dEs_p*(((1/3)*M*m.log(q_f/(3*pf)))+1) ;print(dEa)
    ####
    K=(v*p_f)/Kappa
    E=3*K*(1-2*nu)
    #G=E/(2*(1+nu)) #;print(G)
    p_i=p_f
    q_i=q_f
    p=p+[p_i]
    q=q+[q_i]
    d_evole=d_evole+[dEvol_e]
print(d_evole)
#Dibuja la linea de estado critico
pLEC=np.array(range(1,300,1))
qLEC=M*pLEC
#Dibuja la supericie linea de falla de ROSCOE
pSELR=np.array(range(int(pf),440,1))
qSELR=M*pSELR*(1-(np.log((pSELR)/pf)))
#Dibuja la trayectoria de esfuerzos no drenada
plt.plot(p_1, q_1)
plt.plot(pLEC, qLEC)
plt.plot(pSELR, qSELR)
plt.xlabel("$p' (kPa)$", fontsize = 20)
plt.ylabel("$q (kPa)$", fontsize = 20)
plt.grid(True)
plt.show()