Hi All,
Has anyone (recently) used SAniSandMS successfully in opensees/ openseepy?
The below script is for a single element drained triaxial. It works OK for the older version of Sanisand04 (Manzari Dafalias model), but will not work with SAniSandMS. There is sometimes an error:
SAniSandMS::setResponse -- Unrecognized response option "fabric"
In addition to the above (which could certainly be my modelling error!), if you try the tcl example for SAniSandMS provided in opensees (3.5) – copy and paste it, the output does not resemble the suggested visual output. The cyclic loading is applied properly, but the principle strain generation does not accumulate (ratchet), but cycles at a strain increment of around -4e-3. This is different from the ratcheting displayed in the output of the example and expected performance of the model. This tcl example is located:
https://opensees.github.io/OpenSeesDocu ... andMS.html
So I would be interested to hear if anyone has successfully used SAniSandMS in Openseepy recently?
And, if they have any suggestions ( I am trying very small time steps) …. or even perhaps there is something wrong with the model at the moment/ function parameter list (fabric deprecated?), as in the tcl example, as far I can see the output does not seem to produce the results shown in the example.
Many thanks in advance!
The drained TXL script for Manzari Dafalias (working) / SAniSandMS (strange results) is below (NB long running time due to small dT):
import openseespy.opensees as op
import time
import numpy as np
import matplotlib.pyplot as plt
# Test Specific parameters
pConf = -300.0
devDisp = -0.3
perm = 1.0e-10
vR = 0.8
damp = 0.1
omega1 = 0.0157
omega2 = 64.123
a1 = 2.0 * damp / (omega1 + omega2)
a0 = a1 * omega1 * omega2
op.wipe()
# Create a 3D model with 4 Degrees of Freedom
op.model('BasicBuilder', '-ndm', 3, '-ndf', 4)
# Create nodes
op.node(1, 1.0, 0.0, 0.0)
op.node(2, 1.0, 1.0, 0.0)
op.node(3, 0.0, 1.0, 0.0)
op.node(4, 0.0, 0.0, 0.0)
op.node(5, 1.0, 0.0, 1.0)
op.node(6, 1.0, 1.0, 1.0)
op.node(7, 0.0, 1.0, 1.0)
op.node(8, 0.0, 0.0, 1.0)
# Create Fixities (dof 4 = 1 = drained)
op.fix (1, 0, 1, 1, 1)
op.fix (2, 0, 0, 1, 1)
op.fix (3, 1, 0, 1, 1)
op.fix (4, 1, 1, 1, 1)
op.fix (5, 0, 1, 0, 1)
op.fix (6, 0, 0, 0, 1)
op.fix (7, 1, 0, 0, 1)
op.fix (8, 1, 1, 0, 1)
# Create material
G0 = 110.0 # [Adimensional]
nu = 0.05 # [Adimensional]
e_init = 0.72 # [Adimensional]
Mc = 1.27 # [Adimensional]
c = 0.712 # [Adimensional]
lambda_c = 0.049 # [Adimensional]
e0 = 0.845 # [Adimensional]
ksi = 0.27 # [Adimensional]
P_atm = 101.3 # [kPa]
m = 0.01 # [Adimensional]
h0 = 5.95 # [Adimensional]
ch = 1.01 # [Adimensional]
nb = 2.0 # [Adimensional]
A0 = 1.06 # [Adimensional]
nd = 1.17 # [Adimensional]
z_max = 4 # For SAniSand [Adimensional]
cz = 0 # For SAniSand [Adimensional]
mu0 = 260.0 # For SAniSand [Adimensional]
zeta = 0.0005 # For SAniSand [Adimensional]
beta = 1.0 # For SAniSand [Adimensional]
w1 = 0.5
w2 = 2
#fabric_flag
#flow_flag
Den = 1.584 # [Mg/m^3]
intScheme = 3 # Corresponds to Modified-Euler integration scheme
TanType = 1 # 0: elastic stiffness, 1: continuum elastoplastic stiffness
JacoType = 1 # Not used in explicit methods
TolF = 1.0e-6 # Tolerances, not used in explicit
TolR = 1.0e-6 # Tolerances, not used in explicit
#Reference atmospheric pressure
P_ref = -101.3
# Create material
#nDMaterial SAniSandMS $matTag $G0 $nu $e_init $Mc $c $lambda_c $e0 $ksi $P_atm $m $h0 $ch $nb $A0 $nd $zeta $mu0 $beta $Den $fabric_flag $flow_flag $intScheme $TanType $JacoType $TolF $TolR
#nDMaterial SAniSandMS 1 $G0 $nu $e_init $Mc $c $lambda_c $e0 $ksi $P_atm $m $h0 $ch $nb $A0 $nd $zeta $mu0 $beta $Den $intScheme $TanType $JacoType $TolF $TolR
op.nDMaterial('SAniSandMS', 1, G0, nu, e_init, Mc, c, lambda_c, e0, ksi, P_atm, m, h0, ch, nb, A0, nd, zeta, mu0, beta, Den, intScheme, TanType, JacoType, TolF, TolR)
op.type = "RK"
# nDmaterial ManzariDafalias $matTag $G0 $nu $e_init $Mc $c $lambda_c $e0 $ksi $P_atm $m $h0 $ch $nb $A0 $nd $z_max $cz $Den
#op.nDMaterial('ManzariDafalias', 1, 125, 0.05, vR, 1.25, 0.712, 0.019, 0.934, 0.7, 100, 0.01, 7.05, 0.968, 1.1, 0.704, 3.5, 4, 600, 1.42)
# Create element
op.element('SSPbrickUP', 1, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2.2e6, 1.0, perm, perm, perm, vR, 1.5e-9)
# Create recorders
op.recorder('Node', '-file', 'disp.txt', '-time', '-nodeRange', 1, 8, '-dof', 1, 2, 3, 'disp')
op.recorder('Node', '-file', 'press.txt', '-time', '-nodeRange', 1, 8, '-dof', 4, 'vel')
op.recorder('Element','-file', 'Strain.txt','-time','-ele', 1, 'strain')
op.recorder('Element','-file', 'Stress.txt','-time','-ele', 1, 'stress')
op.recorder('Element','-file', 'Alpha.txt','-time','-ele', 1, 'alpha')
op.recorder('Element','-file', 'Fabric.txt','-time','-ele', 1, 'fabric')
# Create analysis
op.constraints('Transformation')
op.test('NormDispIncr', 1.0e-5, 20, 0)
op.algorithm('Newton')
op.numberer('RCM')
op.system('BandGeneral')
op.integrator('Newmark', 0.5, 0.25)
op.analysis('Transient')
# Apply confinement pressure
pNode = pConf / 4.0
op.timeSeries('Path', 1, '-values', 0, 1, 1, '-time', 0.0, 1000.0, 10000000000.)
op.pattern("Plain", 1, 1, '-factor',1.0)
op.load (1, pNode, 0.0, 0.0, 0.0),
op.load (2, pNode, pNode, 0.0, 0.0),
op.load (3, 0.0, pNode, 0.0, 0.0),
op.load (4, 0.0, 0.0, 0.0, 0.0),
op.load (5, pNode, 0.0, pNode, 0.0),
op.load (6, pNode, pNode, pNode, 0.0),
op.load (7, 0.0, pNode, pNode, 0.0),
op.load (8, 0.0, 0.0, pNode, 0.0)
# attempt at psuedo static loading
op.analyze(10000,0.1)
# Read vertical displacement of top plane
vertDisp = op.nodeDisp(5, 3)
# Apply deviatoric strain
lValues = [1, 1 + devDisp/vertDisp, 1 + devDisp/vertDisp]
# Define the time series
op.timeSeries('Path', 2, '-values', *lValues, '-time', 1000.0, 2000.0, 10000000, '-factor', 1.0)
# Create the loading pattern
op.pattern('Plain', 2, 2, '-factor', 1)
op.sp(5, 3, vertDisp)
op.sp(6, 3, vertDisp)
op.sp(7, 3, vertDisp)
op.sp(8, 3, vertDisp)
# Set number and length of (pseudo)time steps
dT = 0.001
numStep = 1000000
# Initialize the remaining steps and success flag
remStep = numStep
success = 0
def subStepAnalyze(dT, subStep):
if subStep > 10:
return -10
for i in range(1, 3):
print(f"Try dT = {dT}")
success = op.analyze(1, dT)
if success != 0:
success = subStepAnalyze(dT/2.0, subStep + 1)
if success == -10:
print("Did not converge.")
return success
else:
if i == 1:
print(f"Substep {subStep} : Left side converged with dT = {dT}")
else:
print(f"Substep {subStep} : Right side converged with dT = {dT}")
return success
import time
print("Start analysis")
startT = time.time()
while success != -10:
subStep = 0
success = op.analyze(remStep, dT)
if success == 0:
print("Analysis Finished")
break
else:
curTime = op.getTime()
print(f"Analysis failed at {curTime}. Try substepping.")
success = subStepAnalyze(dT / 2.0, subStep + 1)
curStep = int((curTime - 20000) / dT + 1)
remStep = int(numStep - curStep)
print(f"Current step: {curStep}, Remaining steps: {remStep}")
endT = time.time()
print(f"Loading analysis execution time: {endT - startT} seconds.")
op.wipe()
#%%PostProcessing
import pandas as pd
df_stress = pd.read_csv('stress.txt', sep=" ", header=None)
df_strain = pd.read_csv('strain.txt', sep=" ", header=None)
#df_PP = pd.read_csv('CycPP.txt', sep=" ", header=None)
#stress yy = stress xx
Stress_yy = df_stress.iloc[:, 2].to_numpy()*(-1.0) #compression is positive
Principle_Stress = df_stress.iloc[:, 3].to_numpy()*-1.0
Principle_Strain = df_strain.iloc[:, 3].to_numpy()*-1.0
fig = plt.figure(figsize=(10,18))
ax0 = fig.add_subplot(211)
ax0.plot(Principle_Strain, Principle_Stress, label='stress-strain', linewidth=1.2, color='red')
ax0.set_xlabel("Principle strain,%", fontsize=24)
ax0.set_ylabel("Principle stress, kPa", fontsize=24)
ax0.legend(fontsize=24)
ax1 = fig.add_subplot(212)
ax1.plot(Stress_yy, Principle_Stress, label='Stress Path', linewidth=1.2, color='red')
ax1.set_xlabel("Stress yy = xx, kPa", fontsize=24)
ax1.set_ylabel("Principle stress, kPa", fontsize=24)
ax1.legend(fontsize=24)
plt.show()
plt.close()
SanisandMS - has anyone used recently?
Moderators: silvia, selimgunay, Moderators