14.4.3. PM4Sand model undrained cyclic simple shear elementΒΆ

  1. The original model is from 2D Undrained Cyclic Direct Simple Shear Test Using One Element at University of Washington, Department of Civil and Environmental Eng by Geotechnical Eng Group L. Chen, P. Arduino - Feb 2018.

  2. The Python code is converted by Steve Xu from University of Texas at Austin (zhongzexu@utexas.edu), and shown below, which can be downloaded here.

  1# -*- coding: utf-8 -*-
  2"""
  3Created on Sat Apr  9 17:00:41 2022
  4
  5@author: Zhongze Xu, The University of Texas at Austin
  6
  7Openseespy code to run plane-strain stress-controlled undrained cyclic simple shear element
  8to calibrate pm4sand model
  9Basic Units are m, kN and s unless otherwise specified	
 10
 11original opensees code is from:
 12
 132D Undrained Cyclic Direct Simple Shear Test Using One Element        
 14University of Washington, Department of Civil and Environmental Eng   
 15Geotechnical Eng Group, L. Chen, P. Arduino - Feb 2018                
 16Basic Units are m, kN and s unless otherwise specified				
 17
 18"""
 19# from IPython import get_ipython;   
 20# get_ipython().run_line_magic('reset', '-sf')
 21
 22from datetime import datetime
 23import openseespy.opensees as op
 24import numpy as np
 25import matplotlib.pyplot as plt
 26plt.rcParams["font.family"] = "Times New Roman"
 27
 28#==============================================================================
 29#Input Variables
 30'''
 31nDMaterial('PM4Sand', matTag, Dr, G0, hpo, rho, P_atm, h0, e_max, e_min, 
 32           nb, nd, Ado, z_max, c_z, c_e, phi_cv, nu, g_degr, c_dr, c_kaf, 
 33           Q, R, m_par, F_sed, p_sed)
 34'''
 35atm = -101.325
 36sig_v0 = 2.0* atm #initial vertical stress
 37CSR = 0.2 #cyclic stress ratio
 38Cycle_max = 5 #maximxum number of cycles
 39strain_in = 5.0e-6 #strain increment
 40K0 = 0.5
 41nu = K0/(1+K0) #poisson's ratio
 42devDisp = 0.03 #cutoff shear strain
 43perm = 1e-9 #permeability
 44#==============================================================================
 45
 46#primary parameters
 47Dr = 0.5
 48G0 = 476.0
 49hpo = 0.53 #Contraction rate parameter
 50rho = 1.42 #mass density, KN/m3
 51
 52#secondary parameters
 53P_atm = 101.325
 54# all initial stress dependant parameters have negative default values 
 55# and will be calculated during initialization
 56h0 = -1.0  #Variable that adjusts the ratio of plastic modulus to elastic modulus
 57e_max = 0.8
 58e_min = 0.5
 59e_ini = e_max - (e_max - e_min)*Dr #initial void ratio
 60nb = 0.5 #Bounding surface parameter, nb>=0
 61nd = 0.1 #Dilatancy surface parameter, nd>=0
 62Ado = -1.0
 63#Dilatancy parameter, will be computed at the time of initialization if input value is negative
 64z_max = -1.0 #Fabric-dilatancy tensor parameter
 65c_z = 250.0 #Fabric-dilatancy tensor parameter
 66c_e = -1.0 #Fabric-dilatancy tensor parameter
 67phi_cv = 26.0 #Critical state effective friction angle
 68g_degr = 2.0 #Variable that adjusts degradation of elastic modulus with accumulation of fabric
 69c_dr = -1.0 #Variable that controls the rotated dilatancy surface
 70c_kaf = -1.0 # Variable that controls the effect that sustained static shear stresses have on plastic modulus
 71Q = 10.0 #Critical state line parameter
 72R = 1.5 #Critical state line parameter
 73m_par = 0.01#Yield surface constant
 74F_sed = -1.0#Variable that controls the minimum value the reduction factor of the elastic moduli can get during reconsolidation
 75p_sed = -1.0#Mean effective stress up to which reconsolidation strains are enhanced
 76
 77#%%
 78#Rayleigh Damping Parameters
 79'''
 80rayleigh(alphaM, betaK, betaKinit, betaKcomm)
 81'''
 82damp = 0.02
 83omega1 = 0.2
 84omega2 = 20.0
 85a1 = 2.0*damp/(omega1+omega2) #a1 is alphaM
 86a0 = a1*omega1*omega2 #a0 is betaK
 87#%%
 88#create model
 89#Remove the existing model, important!!! 
 90op.wipe()
 91
 92# set modelbuilder
 93op.model('basic', '-ndm', 2, '-ndf', 3)
 94
 95#model nodes
 96x1 = 0.0
 97y1 = 0.0
 98
 99x2 = 1.0
100y2 = 0.0
101
102x3 = 1.0
103y3 = 1.0
104
105x4 = 0.0
106y4 = 1.0
107
108#create nodes
109
110op.node(1, x1, y1)
111op.node(2, x2, y2)
112op.node(3, x3, y3)
113op.node(4, x4, y4)
114
115#boundary conditions
116op.fix(1, 1, 1, 1)
117op.fix(2, 1, 1, 1)
118op.fix(3, 0, 0, 1)
119op.fix(4, 0, 0, 1)
120op.equalDOF(3,4,1,2) #make node 3 and 4 equal displacement at degrees 1 & 2
121
122#material
123#==================================================================
124#nDMaterial('PM4Sand', matTag, D_r, G_o, h_po, Den, P_atm, h_o, e_max, 
125#e_min, n_b, n_d, A_do, z_max, c_z, c_e, phi_cv, nu, g_degr, c_dr, c_kaf, 
126#Q_bolt, R_bolt, m_par, F_sed, p_sed)
127#==================================================================
128op.nDMaterial('PM4Sand', 1, Dr, G0, hpo, rho, P_atm, h0, e_max, e_min, 
129           nb, nd, Ado, z_max, c_z, c_e, phi_cv, nu, g_degr, c_dr, c_kaf, 
130           Q, R, m_par, F_sed, p_sed)
131
132#element
133op.element('SSPquadUP',1, 1,2,3,4, 1, 1.0, 2.2e6, 1.0, perm, perm, e_ini, 1.0e-5)
134
135#create recorders
136op.recorder('Node','-file', 'Cycdisp.txt','-time', '-node',1,2,3,4,'-dof', 1, 2, 'disp')
137op.recorder('Node','-file', 'CycPP.txt','-time', '-node',1,2,3,4,'-dof', 3, 'vel')
138op.recorder('Element','-file', 'Cycstress.txt','-time','-ele', 1, 'stress')
139op.recorder('Element','-file', 'Cycstrain.txt','-time','-ele', 1, 'strain')
140#%%
141#Analysis officially starts here
142op.constraints('Transformation')
143op.test('NormDispIncr', 1.0e-5, 35, 1)
144op.algorithm('Newton')
145op.numberer('RCM')
146op.system('FullGeneral')
147op.integrator('Newmark', 5.0/6.0, 4.0/9.0)
148op.rayleigh(a1, a0, 0.0, 0.0) #modification
149op.analysis('Transient')
150
151#%%apply consolidation pressure
152pNode = sig_v0/2.0 #put confining pressure evenly on two nodes
153
154# create a plain load pattern with time series 1
155op.timeSeries('Path', 1, '-values', 0, 1, 1, '-time', 0.0, 100.0, 1.0e10)
156op.pattern("Plain", 1, 1, '-factor',1.0)
157op.load(3, 0.0, pNode, 0.0) #apply vertical pressure at y direction
158op.load(4, 0.0, pNode, 0.0)
159op.updateMaterialStage('-material', 1, '-stage', 0)
160op.analyze(100,1.0)
161vDisp = op.nodeDisp(3,2)
162b = op.eleResponse(1, 'stress') #b = [sigmaxx, sigmayy, sigmaxy]
163print('shear stress is',b[2])
164op.timeSeries('Path', 2, '-values', 1.0, 1.0, 1.0, '-time', 100.0, 80000.0, 1.0e10, '-factor', 1.0)
165op.pattern('Plain', 2, 2,'-factor',1.0)
166op.sp(3, 2, vDisp)
167op.sp(4, 2, vDisp)
168
169#Close Drainage
170for i in range(4):
171    op.remove('sp', i+1, 3)
172    print('Node ID', i+1)
173
174
175op.analyze(25,1.0)
176b = op.eleResponse(1, 'stress') #b = [sigmaxx, sigmayy, sigmaxy]
177print('shear stress is',b[2])
178print('Drainage is closed')
179
180op.updateMaterialStage('-material', 1, '-stage', 1)
181'''
182Note:
183The program will use the default value of a secondary parameter if 
184a negative input is assigned to that parameter, e.g. Ado = -1. 
185However, FirstCall is mandatory when switching from elastic to elastoplastic 
186if negative inputs are assigned to stress-dependent secondary parameters, 
187e.g. Ado and zmax. 
188'''
189
190#setParameter('-value', 0, '-ele', elementTag, 'FirstCall', matTag)
191op.setParameter('-val', 0, '-ele', 1, 'FirstCall', '1')
192
193op.analyze(25,1.0)
194b = op.eleResponse(1, 'stress') #b = [sigmaxx, sigmayy, sigmaxy]
195print('shear stress is',b[2])
196print('finished update fixties')
197# update Poisson's ratio for analysis
198#setParameter -value 0.3 -ele 1 poissonRatio 1
199op.setParameter('-val', 0.3, '-ele',1, 'poissonRatio', '1')
200
201
202controlDisp = 1.1 * devDisp
203numCycle = 0.25
204print('Current Number of Cycle:', numCycle)
205
206start = datetime.now()
207while (numCycle <= Cycle_max):
208    #impose 1/4 cycle: zero stress to positve max stress  
209    hDisp = op.nodeDisp(3,1)
210    cur_time = op.getTime()
211    steps = controlDisp/strain_in
212    time_change = cur_time + steps
213    op.timeSeries('Path', 3,'-values', hDisp, controlDisp, controlDisp, '-time', cur_time, time_change, 1.0e10, '-factor', 1.0)
214    op.pattern('Plain', 3, 3, '-fact', 1.0)
215    op.sp(3, 1, 1.0)
216    b = op.eleResponse(1, 'stress') #b = [sigmaxx, sigmayy, sigmaxy]
217    print('shear stress is',b[2])
218    while b[2] <= CSR*sig_v0*(-1.0): #b[2] is the shear stress, sigmaxy
219        op.analyze(1, 1.0)
220        b = op.eleResponse(1, 'stress')
221        hDisp = op.nodeDisp(3,1)
222        if hDisp >= devDisp:
223            print('loading break')
224            break
225    numCycle = numCycle + 0.25
226    hDisp = op.nodeDisp(3,1)
227    cur_time = op.getTime()
228    op.remove('loadPattern', 3)
229    op.remove('timeSeries', 3)
230    op.remove('sp', 3, 1)
231    #impose 1/2 cycle: Postive max stress to negative max stress
232    steps = (controlDisp+hDisp)/strain_in
233    time_change = cur_time + steps
234    op.timeSeries('Path', 3,'-values', hDisp, -1.0*controlDisp, -1.0*controlDisp, '-time', cur_time, time_change, 1.0e10, '-factor', 1.0)
235    op.pattern('Plain', 3, 3)
236    op.sp(3, 1, 1.0)
237    while b[2] > CSR*sig_v0:
238        op.analyze(1, 1.0)
239        b = op.eleResponse(1, 'stress')
240        print('shear stress is',b[2])
241        hDisp = op.nodeDisp(3,1)
242        if hDisp <= -1.0*devDisp:
243            print('unloading break')
244            break
245    numCycle = numCycle + 0.5
246    hDisp = op.nodeDisp(3,1)
247    cur_time = op.getTime()
248    op.remove('loadPattern', 3)
249    op.remove('timeSeries', 3)
250    op.remove('sp', 3, 1)
251    #impose 1/4 cycle
252    steps = (controlDisp+hDisp)/strain_in
253    op.timeSeries('Path', 3,'-values', hDisp, controlDisp, controlDisp, '-time', cur_time, time_change, 1.0e10, '-factor', 1.0)
254    op.pattern('Plain', 3, 3, '-fact', 1.0)
255    op.sp(3, 1, 1.0)
256    while b[2] <= 0.0: #b[2] is the shear stress, sigmaxy
257        op.analyze(1, 1.0)
258        b = op.eleResponse(1, 'stress')
259        print('shear stress is',b[2])
260        hDisp = op.nodeDisp(3,1)
261        if hDisp >= devDisp:
262            print('reloading break')
263            break
264    numCycle = numCycle + 0.25
265    print('Current Number of Cycle:', numCycle)
266    op.remove('loadPattern', 3)
267    op.remove('timeSeries', 3)
268    #op.remove('sp', 3, 1)
269
270op.wipe()
271print('Analysis is done!')
272end = datetime.now()
273run_time = end-start
274print('Computation time is' , run_time)
275#%%PostProcessing
276import pandas as pd
277df_stress = pd.read_csv('Cycstress.txt', sep=" ", header=None)
278df_strain = pd.read_csv('Cycstrain.txt', sep=" ", header=None)
279#df_PP = pd.read_csv('CycPP.txt', sep=" ", header=None)
280
281Stress_V = df_stress.iloc[:, 1].to_numpy()*(-1.0) #compression is positive 
282Shear_Stress = df_stress.iloc[:, 3].to_numpy()
283Shear_Strain = df_strain.iloc[:, 3].to_numpy()*100.0
284
285fig = plt.figure(figsize=(10,18))
286ax0 = fig.add_subplot(211)
287ax0.plot(Shear_Strain, Shear_Stress, label='stress-strain', linewidth=0.8)
288ax0.set_xlabel("Shear Strain,%", fontsize=16)
289ax0.set_ylabel("Shear Stress, kPa", fontsize=16)
290ax0.legend(fontsize=16)
291
292ax1 = fig.add_subplot(212)
293ax1.plot(Stress_V, Shear_Stress, label='Stress Path', linewidth=0.8)
294ax1.set_xlabel("Vertical Stress, kPa", fontsize=16)
295ax1.set_ylabel("Shear Stress, kPa", fontsize=16)
296ax1.legend(fontsize=16)
297plt.show()
298plt.close()