################################################
# JOINT_MECA_ENDO: Shear test #
################################################
from math import*
Import numpy as nP
import scipy.optimize as opt
import matplotlib.pyplot as plt
###############################################################
SELECTION DES DONNEES FROM ENTREE:
###############################################################
# Fashion = ['Ilaria', 'CIH']
fashion = 'CIH'
if mode == 'Ilaria':
# Friction coefficient
mu = 1. # Excerpt from ssnap162a line 35
Bn = 137622240.236444 # Excerpt from lcejdm line 166
Bt = 32281760.05546218 # Excerpt from lcejdm line 167
# Energy dissipated during damage:
D1 = 1. # Excerpt from lcejdm line 174
if mode == 'CIH':
# Cohesion
C = 1.4e6 # Email extract from CIH on 8/02
#C = 1.e5 # Excerpt from Ilaria's test
# Maximum shear stress
T = 1.1e6 # Email extract from CIH on 8/02
#T = 0.9*C # Excerpt from Ilaria's test
# Friction coefficient
mu = 0.95 # Email extract from CIH on 8/02
# Normal stiffness:
Kn = 1.e12 # Excerpt from ssnap162a line 30
# Tangential stiffness:
Set = 2. *Kn # Excerpt from ssnp162a line 32
# Residual cohesion
c_bar = 0.6
# Constant compressive stress:
sigma_n_comp = 0.4e6
# Peak shape parameters
m1 = 3. # Excerpt from lcejdm line 172
m2 = 0.5 # Excerpt from lcejdm line 173
# No integration (alpha from 0 to 1 per step step)
step=0.0001
# Max damage (between 0 and 1)
endo_max = 0.90
###############################################################
# Definition of S and R functions as a function of alpha:
###############################################################
# Definition of the S function (alpha):
def S (alpha):
return ((1-alpha) **m1)/(alpha**m2)
# Definition of the function derived from S (alpha):
def Sp (alpha):
return ((m2-m1) *alpha-m2)* ((1-alpha) **(m1-1))/(alpha** (m2+1))
# Definition of the second derivative of S (alpha):
def S2p (alpha):
return ((m1-m2) *(m1-m2-1.)*alpha**2+2.*m2*(m1-m2-1.)*alpha+
m2*(m2+1.))* (1.-alpha**(m1-2)) /alpha** (m2+2)
# Definition of the R function (alpha):
def R (alpha):
return (2). *(Sp (alpha))**2)/(S2p (alpha))) -S (alpha)
#############################################################
# Finding the max of R (alpha):
#############################################################
def findMaxr ():
val_max = opt.fmin (lambda alpha: -R (alpha), 1.e-10)
return R (val_max [0])
#############################################################
# Calculating Bn, Bt, C, and T:
#############################################################
def calcBN ():
Bn = (-Sp (alpha_max) *T**2)/(2. *D1*S (alpha_max) **2)
Return Bn
def calcBT ():
Bt = (-Sp (alpha_max) *(C**2-mu**2*t**2))/(2.*D1*S (alpha_max)**2)
Return Bt
def calCC ():
C = sqrt (2. *D1* (Bn*(mu**2) +Bt))) *S (alpha_max)/(sqrt (-Sp (alpha_max))))
Return C
def calculation ():
T = sqrt (2. *D1*Bn) *S (alpha_max) /sqrt (-Sp (alpha_max))
Return T
#############################################################
# Calculation of delta and sigma as a function of alpha:
#############################################################
def join_meca_endo (Kn, Kt, Bn, Bt, Bn, Bt, m1, m2, mu, c_bar, D1, alpha):
delta_t = (sqrt (2. *D1/ (- (Mu**2*BN+Bt)*Sp (alpha)))) *(Kt+ (Mu**2*BN+Bt)*S (alpha)) /Kt+
(mu*sigma_n_comp+c_bar) /Kt)
sigma_t = sqrt (2. *D1* (mu**2*bn+bt)) *S (alpha) /sqrt (-Sp (alpha)) +mu*sigma_n_comp+c_bar
delta_n = -sigma_n_comp/kn+mu*sqrt (2.*D1/ (- (Mu**2*BN+Bt) *Sp (alpha)))
delta_n_trac = sqrt (2*D1/ (-Bn*Sp (alpha)))) *(Kn+Bn*S (alpha)) /Kn+C_bar/ (Mu*kn)
sigma_n = S (alpha) *sqrt (2.*D1*Bn) /sqrt (-Sp (alpha)) +c_bar/mu
return delta_t, sigma_t, delta_n, delta_n_trac, sigma_n
#############################################################
# Printing the curves:
#############################################################
def plotCurve (x, y, title, label_x, label_x, label_y, output_folder_path, filename):
plt.figure (title)
plt.plot (x, y, label=title)
plt.xlabel (label_x)
plt.ylabel (label_y)
plt.suptitle (title)
plt.legend ()
plt.savefig (output_folder_path+filename, dpi=300)
plt.close ()
if __name__ == '__main__':
# Alpha_max calculation:
alpha_max = (-m2+sqrt (m1*m2/ (m1-m2+1))))/(m1-m2)
# Search for as much R as possible:
R_max = findMaxR ()
# Parameter calculations based on input data:
if mode == 'Ilaria':
C = calCC ()
T = calCT ()
if mode == 'CIH':
D1 = max ((max (T**2/ (2.*Kn), C**2/ (2.*Kt)) *(-Sp (alpha_max)/(S (alpha_max)**2)) *R_max*10.) ,1.) ,1.)
Bn = calcBN ()
Bt = CalcBt ()
# D_min calculation:
D_min = max (T**2/ (2.*Kn), C**2/ (2.*Kt)) *(-Sp (alpha_max)/(S (alpha_max)**2)) *R_max
# Data printing:
Print '\n'
print '===== VALEUR DES PARAMETRES ===== '
print '---------------------------------------------'
Print 'Kn:', Kn, '\ t\ t\ t Kt:', Kt
print 'C:', C, '\ t\ t\ t:', T
print 'mu: ', mu,'\ t\ t\ t c_bar: ', c_bar
print 'Bn: ', Bn,'\ t\ t Bt: ', Bt
print 'alpha_max: ', alpha_max,'\ t R_max: ', R_max:', R_max
print 'D1: ', D1,'\ t\ t sigma_n_comp: ', sigma_n_comp
print 'm1: ', m1,'\ t\ t\ t m2: ', m2
# Parameter validation tests:
Print '\n'
print '===== TEST BY VALIDATION DES PARAMETRES ===== '
print '---------------------------------------------'
print 'Test D1 > D_min:\ t\ t\ t', D1>=d_min
If not D1>D_min:
Print '/! \ D1 = ', D1, 'and D_min =', D_min
print 'Test Mu*t < C < sqrt (2.)*Mu*t:\ t\ t', Mu*t < C < sqrt (2.) *Mu*t
if not mu*T < C < sqrt (2.)*Mute*:
Print '/! \ Mu*t = ', Mu*t, 'C =', C, 'sqrt (2.) *mu*T = ', sqrt (2.) *Mu*t
print 'Test Kn > BN*R_Max:\ t\ t\ t', Kn>BN*R_Max
If not KN>BN*R_Max:
Print '/! \ Kn = ', Kn, 'BN*R_max =', BN*R_max = ', BN*R_max
print 'Test Kt > BT*R_Max:\ t\ t\ t', kt>BT*R_Max
If not KT>BT*R_Max:
Print '/! \ Kt = ', Kt, 'BT*R_Max =', BT*R_Max
Print 'Test Bn > 0\ t\ t\ t\ t', Bn>0.
If not Bn>0. :
Print '/! \ Bn = ', Bn
print 'Test Bt > 0\ t\ t\ t\ t', Bt>0.
If not Bt>0. :
Print '/! \ Bt = ', Bt
print '---------------------------------------------'
Print '\n'
# Delta and sigma calculations:
list_delta_t = [:ref:`0. <0.>`]
list_sigma_t = [:ref:`0. <0.>`]
list_delta_n = [:ref:`0. <0.>`]
list_delta_n_trac = [:ref:`0. <0.>`]
list_sigma_n = [:ref:`0. <0.>`]
for alpha in np.arange (0. +step, endo_max+step, step):
delta_t, sigma_t, delta_n, delta_n_trac, sigma_n = joint_meca_endo (Kn, Kt, Bn, Bt, m1, m2, m2, m2, mu, c_bar, D1, alpha)
list_delta_t.append (delta_t*1e3)
list_sigma_t.append (sigma_t/1.e6)
list_delta_n.append (delta_n*1e3)
list_delta_n_trac.append (delta_n_trac*1e3)
list_sigma_n.append (sigma_n/1.e6)