Thank you Prof. Scott
Search found 10 matches
- Fri Feb 10, 2023 12:42 pm
- Forum: OpenSeesPy
- Topic: Sine Wave load on nodes
- Replies: 2
- Views: 2492
- Mon Feb 06, 2023 2:20 pm
- Forum: OpenSeesPy
- Topic: Sine Wave load on nodes
- Replies: 2
- Views: 2492
Sine Wave load on nodes
Hi all,
I am tranferring Opensees Tcl to Openseepy but encountered a problem here.
The code is to model a 9-4 quad u-p element under cyclic load (sine) with PDMY03 model.
The original code can be found https://opensees.berkeley.edu/wiki/inde ... mentdriver
and I am working on the "undrained cyclic" case
the part of tcl i am having issue with is:
so my python script correspondingly looks like this
I did equalDOF of nodes 3, 4, 7 in DOF 1 and 2, nodes 3 and 4 also equalDOF at degree 3 (fluid pressure)
However, the error i got is:
Any help is highly appreciated,
Thank you.
Regards,
Steve
I am tranferring Opensees Tcl to Openseepy but encountered a problem here.
The code is to model a 9-4 quad u-p element under cyclic load (sine) with PDMY03 model.
The original code can be found https://opensees.berkeley.edu/wiki/inde ... mentdriver
and I am working on the "undrained cyclic" case
the part of tcl i am having issue with is:
Code: Select all
model basic -ndm 2 -ndf 3
pattern Plain 4 "Sine 0 [expr $period*1000] $period" {
load 3 [expr $csr*$vertPress*$xsize*$thickness*0.25] 0 0
load 4 [expr $csr*$vertPress*$xsize*$thickness*0.25] 0 0
}
model basic -ndm 2 -ndf 2
pattern Plain 5 "Sine 0 [expr $period*1000] $period" {
load 7 [expr $csr*$vertPress*$xsize*$thickness*0.5] 0
}
puts "$Analysis_case"
while {[expr abs([nodeDisp 7 1]/$ysize)]<$target_shear_strain} {
# while {[getTime]<0.25}
analyze 1 $deltaT [expr $deltaT/100] $deltaT 20
puts "time = [getTime] sec"
}
Code: Select all
op.constraints('Plain')
op.test('NormDispIncr', 1.0e-5, 50, 0)
op.algorithm('KrylovNewton')
op.numberer('RCM')
op.system('ProfileSPD')
op.integrator('Newmark', 5.0/6.0, 4.0/9.0)
op.rayleigh(Mproprayleigh, 0.0, Kinitproprayleigh, 0.0) #modification
op.analysis('VariableTransient')
op.model('basic', '-ndm', 2, '-ndf', 3)
op.timeSeries('Trig', 3, 0.0, period*1000.0, period, '-factor', 1.0, '-shift', 0.0)
op.pattern('Plain', 3, 3, '-fact', 1.0)
op.sp(4, 1, csr*load_conf)
op.sp(3, 1, csr*load_conf)
op.model('basic', '-ndm', 2, '-ndf', 2)
op.timeSeries('Trig', 4, 0.0, period*1000.0, period, '-factor', 1.0, '-shift', 0.0)
op.pattern('Plain', 4, 4, '-fact', 1.0)
op.sp(7, 1, 2.0*csr*load_conf)
target_shear_strain = 0.03
while abs(op.nodeResponse(3,1,1))/0.1 < target_shear_strain:
op.analyze(1, deltaT, deltaT/100, deltaT, 20)
However, the error i got is:
Code: Select all
WARNING PlainHandler::handle() - non-homogeneos constraint for node 4 homogeneous constraint assumed
WARNING PlainHandler::handle() - non-homogeneos constraint for node 3 homogeneous constraint assumed
WARNING PlainHandler::handle() - non-homogeneos constraint for node 7 homogeneous constraint assumed
WARNING PlainHandler::handle() - constraint at dof 0 already specified for constrained node in MP_Constraint at node 4
WARNING PlainHandler::handle() - constraint at dof 0 already specified for constrained node in MP_Constraint at node 7
Thank you.
Regards,
Steve
- Wed Aug 17, 2022 5:11 pm
- Forum: OpenSeesPy
- Topic: cannot find function updateMaterials
- Replies: 0
- Views: 11051
cannot find function updateMaterials
Dear all,
I am trying to use updateMaterials function to update bulk modulus of UCSD soil models in openseespy but it seems there is no corresponding function called updateMaterials like the one in Opensees (https://opensees.berkeley.edu/OpenSees/ ... l/4246.htm). I cannot find it in the openseespy manual (https://openseespydoc.readthedocs.io/en ... index.html).
The tcl script is:
updateMaterials -material $matTag bulkModulus [expr $G*2/3.]
How do I do this in openseespy?
Any help is greatly appreciated.
Regards,
Steve
I am trying to use updateMaterials function to update bulk modulus of UCSD soil models in openseespy but it seems there is no corresponding function called updateMaterials like the one in Opensees (https://opensees.berkeley.edu/OpenSees/ ... l/4246.htm). I cannot find it in the openseespy manual (https://openseespydoc.readthedocs.io/en ... index.html).
The tcl script is:
updateMaterials -material $matTag bulkModulus [expr $G*2/3.]
How do I do this in openseespy?
Any help is greatly appreciated.
Regards,
Steve
- Sat Jul 30, 2022 7:47 pm
- Forum: OpenSeesPy
- Topic: Use 9_4_quadUP element but python stops working
- Replies: 8
- Views: 3410
Re: Use 9_4_quadUP element but python stops working
Hi all,
Just in case someone is also interested in this topic.
The reason my code didn't work was the I set the wrong boundary condition on two nodes (dof 3, pore pressure).
Just in case someone is also interested in this topic.
The reason my code didn't work was the I set the wrong boundary condition on two nodes (dof 3, pore pressure).
Code: Select all
[/
# -*- coding: utf-8 -*-
"""
Created on Sat Jul 23 11:13:25 2022
@author: Zhongze Xu, the University of Texas at Austin
1D consolidation using fully coupled, u-p formulation.
Cohesive soil. Option of single or double drainage.
Basic units are kN, m, and sec (unless specified)
Details of the model can be found:
https://opensees.berkeley.edu/wiki/index.php/One-dimensional_Consolidation
Original tcl script was written by:
Christopher McGann and Pedro Arduino, University of Washington
"""
#from IPython import get_ipython;
#get_ipython().run_line_magic('reset', '-sf')
from datetime import datetime
import openseespy.opensees as op
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["font.family"] = "Times New Roman"
import sys
from matplotlib import cm
import pandas as pd
#%%
# specify desired drainage type; 1 --> single drainage, 2 --> double drainage
drainage = 1
if drainage == 1:
print('Single Drainage')
elif drainage == 2:
print('Double Drainage')
else:
print('Wrong drainage condition!')
sys.exit(0)
#%%
#create model
#Remove the existing model, important!!!
op.wipe()
#%%
#=========================================
#1. Define Mesh Geometry
#=========================================
# number of elements in horizontal direction
numEleX = 1
# size of elements in horizontal direction (m)
sizeEleX = 1.0
# number of nodes in horizontal direction
numNodeX = numEleX* 2 + 1
# number of elements in vertical direction
numEleY = 10
# size of elements in vertical direction (m)
sizeEleY = 1.0
# number of nodes in vertical direction
numNodeY = 2 * numEleY+1
# total number of nodes
totalNumNode = 3*numNodeY
# total number of elements
totalNumEle = numEleY*numEleX
node_list = []
#%%
#-------------------------------------------------------------------------------------------
# 2. DEFINE CORNER NODES
#-------------------------------------------------------------------------------------------
op.model('basic', '-ndm', 2, '-ndf', 3) #corner nodes record PP
fig = plt.figure(figsize=(10,10))
ax0 = fig.add_subplot(121)
ax0.set_xlabel("x, m", fontsize=16)
ax0.set_ylabel("y, m", fontsize=16)
ax0.set_aspect('equal')
ax0.set_xlim(-1.0, 2.0)
ax0.set_ylim(-1.0,numEleY+1)
for i in range(0,numEleX+2,2):
#print(i)
for j in range(0,2*numEleY+2,2):
#print(j)
nodeNum = i+1 + (j)*numNodeX
node_list.append(nodeNum)
Xcoord = (i)*sizeEleX/2
Ycoord = (j)*sizeEleX/2
#set corner nodes
op.node(nodeNum, Xcoord, Ycoord)
#plot corner nodes
ax0.plot(Xcoord, Ycoord, 'ro', label = str(nodeNum))
print('Finished creating corner soil nodes...')
#%%
#-------------------------------------------------------------------------------------------
# 3. DEFINE BOUNDARY CONDITIONS FOR CORNER NODES
#-------------------------------------------------------------------------------------------
# boundary conditions for base nodes depending upon selected drainage type
if drainage == 1: #single drainage, base deformation fixed, pore pressure free = no drainage
op.fix(1, 1, 1, 0)
op.fix(3, 1, 1, 0)
elif drainage == 2: #double drainages, base deformation fixed, pore pressure fixed, free drainage
op.fix(1, 1, 1, 1)
op.fix(3, 1, 1, 1)
#change the surface node to free-drainage no displacement in both directions
op.fix(totalNumNode, 1 , 0, 1)
op.fix(totalNumNode-2, 1 , 0, 1)
#print('Nodes with fixed X and Y displacement are:')
for i in range(numEleY): #first fix all other element in the x direction
if i != 0:
nodeNum3 = 6*(i+1)-5
nodeNum4 = 6*(i+1)-3
op.fix(nodeNum3, 1, 0, 0)
op.fix(nodeNum4, 1, 0, 0)
print(nodeNum3)
print(nodeNum4)
print('Finished creating boundary conditions for corner nodes...')
#%%
#-------------------------------------------------------------------------------------------
# 4. DEFINE INTERIOR NODES
#-------------------------------------------------------------------------------------------
op.model('basic', '-ndm', 2, '-ndf', 2) #interior nodes do not record PP, 2 dof
for i in range(numEleY):
nodeNum0 = 6*(i+1)-4
nodeNum1 = 6*(i+1)-2
nodeNum2 = 6*(i+1)-1
nodeNum3 = 6*(i+1)
Xcoord0 = 0.5
Xcoord1 = 0.0
Xcoord2 = 0.5
Xcoord3 = 1.0
Ycoord0 = i*1.0
Ycoord1 = i+0.5
Ycoord2 = i+0.5
Ycoord3 = i+0.5
node_list.append(nodeNum0)
node_list.append(nodeNum1)
node_list.append(nodeNum2)
node_list.append(nodeNum3)
print(nodeNum0, nodeNum1, nodeNum2, nodeNum3)
#set interior nodes
op.node(nodeNum0, Xcoord0, Ycoord0)
op.node(nodeNum1, Xcoord1, Ycoord1)
op.node(nodeNum2, Xcoord2, Ycoord2)
op.node(nodeNum3, Xcoord3, Ycoord3)
ax0.plot(Xcoord0, Ycoord0, 'bo', label = str(nodeNum0))
ax0.plot(Xcoord1, Ycoord1, 'bo', label = str(nodeNum1))
ax0.plot(Xcoord2, Ycoord2, 'bo', label = str(nodeNum2))
ax0.plot(Xcoord3, Ycoord3, 'bo', label = str(nodeNum3))
nodeNum_surface3 = 6*(i+1)+2
Xcoord = 0.5
Ycoord = i+1.0
op.node(nodeNum_surface3, Xcoord, Ycoord)
ax0.plot(Xcoord, Ycoord, 'bo', label = str(nodeNum0))
print('Finished creating all soil nodes...')
#%%
#-------------------------------------------------------------------------------------------
# 5. DEFINE BOUNDARY CONDITIONS FOR INTERIOR NODES
#-------------------------------------------------------------------------------------------
op.fix(2, 1, 1)
for i in range(numEleY):
nodeNum0 = 6*(i+1)-4
nodeNum1 = 6*(i+1)-2
nodeNum2 = 6*(i+1)-1
nodeNum3 = 6*(i+1)
if i==0:
op.fix(nodeNum1, 1, 0)
op.fix(nodeNum2, 1, 0)
op.fix(nodeNum3, 1, 0)
#print(nodeNum1,nodeNum2,nodeNum3)
else:
op.fix(nodeNum0, 1, 0)
op.fix(nodeNum1, 1, 0)
op.fix(nodeNum2, 1, 0)
op.fix(nodeNum3, 1, 0)
#print(nodeNum0,nodeNum1,nodeNum2,nodeNum3)
op.fix(nodeNum_surface3, 1, 0)
nodeNum_surface1 = 6*numEleY+3
nodeNum_surface2 = 6*numEleY+1
#make surface node equal displacement at degrees 1 & 2
op.equalDOF(nodeNum_surface1,nodeNum_surface2,1,2)
op.equalDOF(nodeNum_surface1,nodeNum_surface3,1,2)
print('Finished creating all boundary conditions and equalDOF...')
plt.show()
plt.close()
#-------------------------------------------------------------------------------------------
# 6. DEFINE MATERIAL PROPERTIES OF SOIL
#-------------------------------------------------------------------------------------------
# Constutive Model: PressureIndependMultiYield material
# saturated mass density (Mg/m3)
satDensity = 2.3
# fluid mass density (Mg/m3)
H2ODensity = 1.0
# shear modulus (kPa)
shear = 2.5e4
# bulk modulus (kPa)
bulk = 6.2e5
# undrained bulk modulus (kPa)
uBulk = 2.2e5
# vertical permeability, input value in m/s, units are modified to suit constitutive model
mPermV = 5.0e-5
vPerm = mPermV/9.81/H2ODensity
# horizontal permeability, input value in m/s, units are modified to suit constitutive model
mPermH = 5.0e-5
hPerm = mPermH/9.81/H2ODensity
# cohesion (kPa)
cohesion = 45.0
# friction angle (degrees)
phi = 0.0
# peak shear strain
gammaPeak = 0.1
# reference pressure (kPa)
refP = 80.0
# pressure dependency coefficient
pressCoeff = 0.0
# number of yield surfaces
numSurf = 22
#-------------------------------------------------------------------------------------------
# 7. DEFINE SOIL MATERIALS
#-------------------------------------------------------------------------------------------
'''
Soil Model: PressureIndependMultiYield
nDMaterial('PressureIndependMultiYield', matTag, nd, rho, refShearModul, refBulkModul,
cohesi, peakShearStra, frictionAng=0., refPress=100., pressDependCoe=0.,
noYieldSurf=20, *yieldSurf)
'''
matTag = 1
nd = 2 #Number of dimensions, 2 for plane-strain, and 3 for 3D analysis.
rho = satDensity # Saturated Soil Mass Density
op.nDMaterial('PressureIndependMultiYield', matTag, 2, satDensity, shear, bulk,
cohesion, gammaPeak, phi, refP, 0., numSurf)
print('Finished creating all soil materials....')
#-------------------------------------------------------------------------------------------
# 8. DEFINE SOIL ELEMENTS
#-------------------------------------------------------------------------------------------
#element('9_4_QuadUP', eleTag, *eleNodes, thick, matTag, bulk, fmass, hPerm, vPerm, <b1=0, b2=0>)
thick = 1.0
wgtX = 0.0
wgtY = -9.81
'''
i = 0
n0 = 6*i
n1 = n0+1
n2 = n0+2
n3 = n0+3
n4 = n0+4
n5 = n0+5
n6 = n0+6
n7 = n0+7
n8 = n0+8
#op.element('9_4_QuadUP', 1, n0, n1, n2, n3, n4, n5, n6, n7, n8, thick, matTag, bulk, H2ODensity, hPerm, vPerm, wgtX, wgtY)
'''
for i in range(numEleY):
eleTag = i+1
n1 = 6*(i+1)-5
n2 = n1+1
n3 = n1+2
n4 = n1+3
n5 = n1+4
n6 = n1+5
n7 = n1+6
n8 = n1+7
n9 = n1+8
print('Node:', n1, n3, n9, n7, n2, n6, n8, n4, n5)
#Need to check the order of nodes input to the element commend\
#https://opensees.berkeley.edu/wiki/index.php/Four_Node_Quad_u-p_Element
op.element('9_4_QuadUP', eleTag, n1, n3, n9, n7, n2, n6, n8, n4, n5, thick,
matTag, bulk, H2ODensity, hPerm, vPerm, wgtX, wgtY)
# update material stage to zero for elastic behavior
op.updateMaterialStage('-material', matTag, '-stage', 0)
print('Finished creating all soil elements...')
#%%
#-------------------------------------------------------------------------------------------
# 8. APPLY GRAVITY LOADING
#-------------------------------------------------------------------------------------------
gamma = 0.5
beta = 0.25
#Analysis officially starts here
op.constraints('Transformation')
op.test('NormDispIncr', 1.0e-5, 30, 1)
op.algorithm('Newton')
op.numberer('RCM')
op.system('ProfileSPD')
op.integrator('Newmark', gamma, beta)
op.analysis('Transient')
op.analyze(10, 500)
print('Finished with elastic gravity analysis...')
#%%
#-------------------------------------------------------------------------------------------
# 9. UPDATE SOIL MATERIAL STAGE
#-------------------------------------------------------------------------------------------
# update material stage to consider elastoplastic behavior
op.updateMaterialStage('-material', matTag, '-stage', 1)
op.analyze(40, 500)
print('Finished with plastic gravity analysis...')
#%%
#-------------------------------------------------------------------------------------------
# 10. CREATE RECORDERS
#-------------------------------------------------------------------------------------------
op.recorder('Node','-file', 'displacement.out','-time', '-nodeRange',1,totalNumNode,'-dof', 1, 2, 'disp')
op.recorder('Node','-file', 'porepressure.out','-time', '-nodeRange',1,totalNumNode,'-dof', 3, 'vel')
#record stress and strain at the 9th Gaussian point, center of the element
op.recorder('Element','-file','stress.out','-time','-eleRange', 1, numEleY,'material','9','stress')
op.recorder('Element','-file','strain.out','-time','-eleRange', 1, numEleY,'material','9','strain')
print('Finished creating all recorders...')
#%%
#-------------------------------------------------------------------------------------------
# 11. CONSOLIDATION ANALYSIS
#-------------------------------------------------------------------------------------------
setTime = 0.0
op.wipeAnalysis()
overburden = -1000.0 #set vertical load
op.timeSeries('Constant', 1) #tsTag = 1
op.pattern('Plain', 1, 1)
op.load(totalNumNode-1, 0.0,overburden) #put overburden on y direction
print('Consolidation loading created...')
op.constraints('Penalty', 1.0e16, 1.0e16)
op.test('NormDispIncr', 1.0e-4, 30, 1)
op.algorithm('Newton')
op.numberer('RCM')
op.system('ProfileSPD')
op.integrator('Newmark', gamma, beta)
op.rayleigh(0.5, 0.5, 0.0, 0.0)
op.analysis('Transient')
if drainage == 1:
op.analyze(800, 0.25)
else:
op.analyze(2000, 0.25)
print('Finished with consolidation analysis...')
op.wipe()
#%%
#-------------------------------------------------------------------------------------------
# 12. Post-Process
#-------------------------------------------------------------------------------------------
df = pd.read_csv('porepressure.out', delimiter=' ', header=None)
df_ID = df.iloc[:,0]
step_ID = df_ID.to_numpy()
step_ID = step_ID-step_ID[0]
df = df.iloc[: , 1:]
nStep = df.shape[0]
nNode = df.shape[1]
plot_times = nStep
DOC = step_ID/step_ID[-1] #something like degree of consolidation for plot only!!!!
depth_left = np.linspace(-numEleY, 0, numEleY+1)
fig= plt.figure(figsize=(6,6))
ax0 = fig.add_subplot(111)
ax0.set_xlabel('Excess Pore Pressure, kPa', fontsize=20, fontweight='bold')
ax0.set_ylabel('Depth, m', fontsize=20, fontweight='bold')
ax0.set_ylim(-numEleY-0.5,0.5)
ax0.tick_params(direction='in', labelsize=13)
ax0.tick_params(axis='both', which='both',direction='in',labelsize=14, width=1.5, length=3)
for i in range(plot_times):
plot_arr = df.iloc[i,:].to_numpy()
PP_left = plot_arr[0:nNode:6]
DOC_arr = DOC[i]*np.ones(numEleY+1)
if i > 5:
im = ax0.plot(PP_left,depth_left, lw=3, c=cm.turbo(DOC[i]))
if drainage == 1:
ax0.text(0.9, 0.9,'Single Drainage', fontsize = 18,horizontalalignment='right',
verticalalignment='top',transform = ax0.transAxes)
else:
ax0.text(0.9, 0.9,'Double Drainage', fontsize = 18,horizontalalignment='right',
verticalalignment='top',transform = ax0.transAxes)
plt.show()
plt.close()
- Thu Jul 28, 2022 9:14 am
- Forum: OpenSeesPy
- Topic: Use 9_4_quadUP element but python stops working
- Replies: 8
- Views: 3410
- Wed Jul 27, 2022 10:08 pm
- Forum: OpenSeesPy
- Topic: Use 9_4_quadUP element but python stops working
- Replies: 8
- Views: 3410
- Wed Jul 27, 2022 8:18 pm
- Forum: OpenSeesPy
- Topic: Use 9_4_quadUP element but python stops working
- Replies: 8
- Views: 3410
Re: Use 9_4_quadUP element but python stops working
sorry that the original code I posted has some errors but I have update the code.
If there is wrong with the code, error or warning should be given by python but the only thing I can see is the warning that kernel dead.
I tried with Spyder and Jupyter Notebook. Both Macos and Windows 11. All same issue.
Any help is greatly appreciated!
Thanks in advance.
Regards
Steve
If there is wrong with the code, error or warning should be given by python but the only thing I can see is the warning that kernel dead.
I tried with Spyder and Jupyter Notebook. Both Macos and Windows 11. All same issue.
Any help is greatly appreciated!
Thanks in advance.
Regards
Steve
- Mon Jul 25, 2022 12:59 pm
- Forum: OpenSeesPy
- Topic: Use 9_4_quadUP element but python stops working
- Replies: 8
- Views: 3410
Use 9_4_quadUP element but python stops working
I am practicing a 1-D consolidation with 9_4_quadUP elements.
The original tcl code can be found at https://opensees.berkeley.edu/wiki/inde ... solidation
The problem I have is that
op.element('9_4_QuadUP', eleTag, nI, nJ, nK, nL, nM, nN, nP, nQ, nR, thick, matTag, bulk, H2ODensity, hPerm, vPerm, wgtX, wgtY)
this line always leads to the restart of python kernel for some reason I do not understand.
The code runs until the line with element command.
Python version: 3.9.12
Any help is greatly appreciated.
Regards,
Steve
Python Code:
The original tcl code can be found at https://opensees.berkeley.edu/wiki/inde ... solidation
The problem I have is that
op.element('9_4_QuadUP', eleTag, nI, nJ, nK, nL, nM, nN, nP, nQ, nR, thick, matTag, bulk, H2ODensity, hPerm, vPerm, wgtX, wgtY)
this line always leads to the restart of python kernel for some reason I do not understand.
The code runs until the line with element command.
Python version: 3.9.12
Any help is greatly appreciated.
Regards,
Steve
Python Code:
Code: Select all
[
from datetime import datetime
import openseespy.opensees as op
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["font.family"] = "Times New Roman"
import sys
#%%
# specify desired drainage type; 1 --> single drainage, 2 --> double drainage
drainage = 1
if drainage == 1:
print('Single Drainage')
elif drainage == 2:
print('Double Drainage')
else:
print('Wrong drainage condition!')
sys.exit(0)
#%%
#create model
#Remove the existing model, important!!!
op.wipe()
#%%
#=========================================
#1. Define Mesh Geometry
#=========================================
# number of elements in horizontal direction
numEleX = 1
# size of elements in horizontal direction (m)
sizeEleX = 1.0
# number of nodes in horizontal direction
numNodeX = numEleX* 2 + 1
# number of elements in vertical direction
numEleY = 11 # numEleY should be an odd number
# size of elements in vertical direction (m)
sizeEleY = 1.0
# number of nodes in vertical direction
numNodeY = 2 * numEleY+1
# total number of nodes
totalNumNode = 3*numNodeY
# total number of elements
totalNumEle = numEleY*numEleX
node_list = []
#-------------------------------------------------------------------------------------------
# 2. DEFINE CORNER NODES
#-------------------------------------------------------------------------------------------
op.model('basic', '-ndm', 2, '-ndf', 3) #corner nodes record PP
fig = plt.figure(figsize=(5,10))
ax0 = fig.add_subplot(211)
ax0.set_xlabel("x, m", fontsize=16)
ax0.set_ylabel("y, m", fontsize=16)
ax0.set_aspect('equal')
ax0.set_xlim(-1.0, 2.0)
ax0.set_ylim(-1.0,numEleY+1)
for i in range(numEleY):
if (i % 2) == 0:
nodeNum0 = 6*i
nodeNum1 = 6*i+2
nodeNum2 = 6*i+6
nodeNum3 = 6*i+8
node_list.append(nodeNum0)
node_list.append(nodeNum1)
node_list.append(nodeNum2)
node_list.append(nodeNum3)
Xcoord0 = 0.0
Xcoord1 = 1.0
Xcoord2 = 0.0
Xcoord3 = 1.0
Ycoord0 = i*1.0
Ycoord1 = i*1.0
Ycoord2 = i*1.0 + 1.0
Ycoord3 = i*1.0 + 1.0
#set corner nodes
op.node(nodeNum0, Xcoord0, Ycoord0)
op.node(nodeNum1, Xcoord1, Ycoord1)
op.node(nodeNum2, Xcoord2, Ycoord2)
op.node(nodeNum3, Xcoord3, Ycoord3)
#plot corner nodes
ax0.plot(Xcoord0, Ycoord0, 'ro', label = str(nodeNum0))
ax0.plot(Xcoord1, Ycoord1, 'ro', label = str(nodeNum1))
ax0.plot(Xcoord2, Ycoord2, 'ro', label = str(nodeNum2))
ax0.plot(Xcoord3, Ycoord3, 'ro', label = str(nodeNum3))
print('Finished creating corner soil nodes...')
nodeNum_surface1 = nodeNum3 #corner node on the surface
nodeNum_surface2 = nodeNum2
print('Number of nodes are:',nodeNum3+1)
#-------------------------------------------------------------------------------------------
# 3. DEFINE BOUNDARY CONDITIONS FOR CORNER NODES
#-------------------------------------------------------------------------------------------
# boundary conditions for base nodes depending upon selected drainage type
if drainage == 1: #single drainage, base fixed
op.fix(0, 1, 1, 0)
op.fix(2, 1, 1, 0)
elif drainage == 2: #double drainages, base deformation fixed, drainage free
op.fix(0, 1, 1, 1)
op.fix(2, 1, 1, 1)
#change the surface node to free-drainage no displacement in both directions
op.fix(nodeNum_surface1, 1 , 1, 0)
op.fix(nodeNum_surface2, 1 , 1, 0)
#print('Nodes with fixed X and Y displacement are:')
for i in range(numEleY): #first fix all other element in the x direction
if (i % 2) == 0:
nodeNum0 = 6*i
nodeNum1 = 6*i+2
nodeNum2 = 6*i+6
nodeNum3 = 6*i+8
if i != 0:
op.fix(nodeNum0, 1, 0, 0)
op.fix(nodeNum1, 1, 0, 0)
#print(nodeNum0)
#print(nodeNum1)
if i != numEleY-1:
op.fix(nodeNum2, 1, 0, 0)
op.fix(nodeNum3, 1, 0, 0)
#print(nodeNum2)
#print(nodeNum3)
print('Finished creating boundary conditions for corner nodes...')
#-------------------------------------------------------------------------------------------
# 4. DEFINE INTERIOR NODES
#-------------------------------------------------------------------------------------------
op.model('basic', '-ndm', 2, '-ndf', 2) #interior nodes do not record PP, 2 dof
for i in range(numEleY):
nodeNum0 = 6*i+1
nodeNum1 = 6*i+3
nodeNum2 = 6*i+4
nodeNum3 = 6*i+5
Xcoord0 = 0.5
Xcoord1 = 0.0
Xcoord2 = 0.5
Xcoord3 = 1.0
Ycoord0 = i*1.0
Ycoord1 = i+0.5
Ycoord2 = i+0.5
Ycoord3 = i+0.5
node_list.append(nodeNum0)
node_list.append(nodeNum1)
node_list.append(nodeNum2)
node_list.append(nodeNum3)
#set interior nodes
op.node(nodeNum0, Xcoord0, Ycoord0)
op.node(nodeNum1, Xcoord1, Ycoord1)
op.node(nodeNum2, Xcoord2, Ycoord2)
op.node(nodeNum3, Xcoord3, Ycoord3)
ax0.plot(Xcoord0, Ycoord0, 'bo', label = str(nodeNum0))
ax0.plot(Xcoord1, Ycoord1, 'bo', label = str(nodeNum1))
ax0.plot(Xcoord2, Ycoord2, 'bo', label = str(nodeNum2))
ax0.plot(Xcoord3, Ycoord3, 'bo', label = str(nodeNum3))
nodeNum_surface3 = 6*i+7
Xcoord = 0.5
Ycoord = i+1.0
op.node(nodeNum_surface3, Xcoord, Ycoord)
ax0.plot(Xcoord, Ycoord, 'bo', label = str(nodeNum0))
print('Finished creating all soil nodes...')
#-------------------------------------------------------------------------------------------
# 5. DEFINE BOUNDARY CONDITIONS FOR INTERIOR NODES
#-------------------------------------------------------------------------------------------
op.fix(1, 1, 1)
for i in range(numEleY):
nodeNum0 = 6*i+1
nodeNum1 = 6*i+3
nodeNum2 = 6*i+4
nodeNum3 = 6*i+5
if i==0:
op.fix(nodeNum1, 1, 0)
op.fix(nodeNum2, 1, 0)
op.fix(nodeNum3, 1, 0)
#print(nodeNum1,nodeNum2,nodeNum3)
else:
op.fix(nodeNum0, 1, 0)
op.fix(nodeNum1, 1, 0)
op.fix(nodeNum2, 1, 0)
op.fix(nodeNum3, 1, 0)
#print(nodeNum0,nodeNum1,nodeNum2,nodeNum3)
op.fix(nodeNum_surface3, 1, 0)
#make surface node equal displacement at degrees 1 & 2
op.equalDOF(nodeNum_surface1,nodeNum_surface2,1,2)
op.equalDOF(nodeNum_surface1,nodeNum_surface3,1,2)
print('Finished creating all boundary conditions and equalDOF...')
#-------------------------------------------------------------------------------------------
# 6. DEFINE MATERIAL PROPERTIES OF SOIL
#-------------------------------------------------------------------------------------------
# Constutive Model: PressureIndependMultiYield material
# saturated mass density (Mg/m3)
satDensity = 2.3
# fluid mass density (Mg/m3)
H2ODensity = 1.0
# shear modulus (kPa)
shear = 2.5e4
# bulk modulus (kPa)
bulk = 6.2e5
# undrained bulk modulus (kPa)
uBulk = 2.2e5
# vertical permeability, input value in m/s, units are modified to suit constitutive model
mPermV = 5.0e-5
vPerm = mPermV/9.81/H2ODensity
# horizontal permeability, input value in m/s, units are modified to suit constitutive model
mPermH = 5.0e-5
hPerm = mPermH/9.81/H2ODensity
# cohesion (kPa)
cohesion = 45.0
# friction angle (degrees)
phi = 0.0
# peak shear strain
gammaPeak = 0.1
# reference pressure (kPa)
refP = 80.0
# pressure dependency coefficient
pressCoeff = 0.0
# number of yield surfaces
numSurf = 22
#-------------------------------------------------------------------------------------------
# 7. DEFINE SOIL MATERIALS
#-------------------------------------------------------------------------------------------
'''
Soil Model: PressureIndependMultiYield
nDMaterial('PressureIndependMultiYield', matTag, nd, rho, refShearModul, refBulkModul,
cohesi, peakShearStra, frictionAng=0., refPress=100., pressDependCoe=0.,
noYieldSurf=20, *yieldSurf)
'''
matTag = 1
nd = 2 #Number of dimensions, 2 for plane-strain, and 3 for 3D analysis.
rho = satDensity # Saturated Soil Mass Density
op.nDMaterial('PressureIndependMultiYield', matTag, 2, satDensity, shear, bulk,
cohesion, gammaPeak, phi, refP, 0., numSurf)
print('Finished creating all soil materials....')
#-------------------------------------------------------------------------------------------
# 8. DEFINE SOIL ELEMENTS
#-------------------------------------------------------------------------------------------
#element('9_4_QuadUP', eleTag, *eleNodes, thick, matTag, bulk, fmass, hPerm, vPerm, <b1=0, b2=0>)
thick = 1.0
wgtX = 0.0
wgtY = -9.81
i = 0
n0 = 6*i
n1 = n0+1
n2 = n0+2
n3 = n0+3
n4 = n0+4
n5 = n0+5
n6 = n0+6
n7 = n0+7
n8 = n0+8
op.element('9_4_QuadUP', i, n0, n1, n2, n3, n4, n5, n6, n7, n8, thick, matTag, bulk, H2ODensity, hPerm, vPerm, wgtX, wgtY)
- Tue May 03, 2022 7:23 am
- Forum: OpenSeesPy
- Topic: Transfer TCL to Openseespy
- Replies: 3
- Views: 2916
Re: Transfer TCL to Openseespy
Thank you very much Prof. Scott!mhscott wrote: ↑Tue May 03, 2022 6:16 am ops.timeSeries('Path',1,'-time',...) #https://openseespydoc.readthedocs.io/en ... athTs.html
ops.pattern('Plain',2,1)
ops.sp(3,...)
ops.sp(4,...)
- Mon May 02, 2022 8:53 pm
- Forum: OpenSeesPy
- Topic: Transfer TCL to Openseespy
- Replies: 3
- Views: 2916
Transfer TCL to Openseespy
Dear all,
I am transferring the opensees pm4sand example to openseespy.
I am currently stuck on one command:
set ts1 "{Series -time {100 80000 1.0e10} -values {1.0 1.0 1.0} -factor 1}"
eval "pattern Plain 2 $ts1 {
sp 3 2 $vDisp
sp 4 2 $vDisp
}"
Openseespy does not seems to have this Series function or in the timeSeries function, it does not have this ramped load. So what is the equivalent function for this in openseespy. I am quite confused.
Any help or suggestion is greatly appreciated.
Regards,
Steve
I am transferring the opensees pm4sand example to openseespy.
I am currently stuck on one command:
set ts1 "{Series -time {100 80000 1.0e10} -values {1.0 1.0 1.0} -factor 1}"
eval "pattern Plain 2 $ts1 {
sp 3 2 $vDisp
sp 4 2 $vDisp
}"
Openseespy does not seems to have this Series function or in the timeSeries function, it does not have this ramped load. So what is the equivalent function for this in openseespy. I am quite confused.
Any help or suggestion is greatly appreciated.
Regards,
Steve