14.4.3. PM4Sand model undrained cyclic simple shear elementΒΆ
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.
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()