I am posting this here, since the soil modelling forum is locked.
I am trying to conduct a soil column response analysis. The elements have a size of 1m X 1m and the column's width is 2m. The flexible rock halfspace is modelled with a viscous damper, hence the type of loading is a Path Type time series of velocity.
My problem is that when i try to model a column with a height greater than 101m, an error occurs that says: "WARNING - PathSeries::PathSeries() - could not open file Velocity1D.out" which is really peculiar, considering the fact that for a height of 100m and less, the analysis carries on as it should. I have no clue as to what the problem is and how I can overcome it. Would gladly welcome any suggestions.
Thank you in advance
Soil column response
Moderators: silvia, selimgunay, Moderators
Re: Soil column response
that error should not occur .. post script.
Re: Soil column response
Since my intention was to be able to modify the characteristics of the column, almost everything is expressed in terms of parameters and nodes and elements are created through loops, so it may be hard to read. The code is as follows:
wipe all;
source INPUT_1D\\OPENSEES1DINPUT.txt
file mkdir OUTPUT_1D\\OS_temp_1D
#-------------------------------------------------------------------------------------------
# INPUT VARIABLES - PRELIMINARY CALCULATIONS
#-------------------------------------------------------------------------------------------
## Parameters that are sourced
#SOILHEIGHT
#SOILVS
#SOILDENSITY
#SOILPOISSON
#cohesion
#phi
#ROCKVS
#ROCKDENSITY
#RECDT
#RECSTEPS
#Define constants - calculate soil parameters - depth for output
set g 9.81
set pi 3.141592654
set thick 1.0; #Element thickness
set ColWidth 2
set WidthMesh 1
set ColWidthMesh [expr $ColWidth/$WidthMesh]
set HeightMesh 1
set ColHeightMesh [expr $SOILHEIGHT/$HeightMesh]
set pressCoeff 0.0
#########################################################
# #
# DEFINE LAYER PROPERTIES #
# #
#########################################################
set d 0.2
set numLayers $ColHeightMesh
set nu $SOILPOISSON
set rho $SOILDENSITY
set cohesion [expr $cohesion*1.0]
set SumVs 0.0
#Set Vsref to match Vs30=SOILVS
for {set k 1} {$k <= 30} {incr k 1} {
set SumVs [expr $SumVs+$SOILVS*pow(($k),$d/2)]
}
set Vs30 [expr $SumVs/30]
set VsRat [expr $SOILVS/$Vs30]
set SOILVSff [expr $SOILVS*$VsRat]
set profileFile [open OUTPUT_1D\\VsProfile.txt w]
# Vs for each layer (m/s)
for {set k $numLayers} {$k >= 1} {incr k -1} {
set Vs($k) [expr $SOILVSff*pow(($numLayers+1-$k),$d/2)]
puts $profileFile "$Vs($k)"
}
close $profileFile
# soil shear modulus for each layer (kPa)
for {set k 1} {$k <= $numLayers} {incr k 1} {
set G($k) [expr $rho*$Vs($k)*$Vs($k)]
}
# soil elastic modulus for each layer (kPa)
for {set k 1} {$k <= $numLayers} {incr k 1} {
set E($k) [expr 2*$G($k)*(1+$nu)]
}
# soil bulk modulus for each layer (kPa)
for {set k 1} {$k <= $numLayers} {incr k 1} {
set bulk($k) [expr $E($k)/(3*(1-2*$nu))]
}
# reference pressure for each layer (kPa)
for {set k 1} {$k <= $numLayers} {incr k 1} {
set refPress($k) 80.0; #Constant
}
# calculate soil period and frequencies
set SOILPERIOD 0.0
for {set k $numLayers} {$k >= 1} {incr k -1} {
set SOILPERIOD [expr $SOILPERIOD + $HeightMesh/$Vs($k)]
}
set SOILPERIOD [expr 4*$SOILPERIOD]
set omega1 [expr 2*$pi / $SOILPERIOD]
set omega2 [expr 5*$omega1]
#-------------------------------------------------------------------------------------------
# 1. DEFINE MODEL BUILDER
#-------------------------------------------------------------------------------------------
model basic -ndm 2 -ndf 2; # 2 spacial dimensions,2 DOF's per node
#----------------------------------------------------------------------------------------------------
# 2. DEFINE SOIL NODES
#----------------------------------------------------------------------------------------------------
set ntag 0
for {set iy 0} {$iy<=$ColHeightMesh} {incr iy 1} {
set nodey [expr $iy*$HeightMesh]
for {set ix 0} {$ix <= $ColWidthMesh} {incr ix 1} {
set nodex [expr $ix*$WidthMesh]
set ntag [expr $ntag+1]
node $ntag $nodex $nodey
}
}
puts "Finished creating soil nodes..."
set NODESNUM $ntag
#-----------------------------------------------------------------------------------------------------------
# 3. DEFINE DASHPOT NODES
#-----------------------------------------------------------------------------------------------------------
node 4000 0.0 0.0
node 4001 0.0 0.0
#-----------------------------------------------------------------------------------------------------------
# 4. DEFINE FIXITY AND EQUAL DOF
#-----------------------------------------------------------------------------------------------------------
# Define fixity of dashpot nodes
fix 4000 1 1
fix 4001 0 1
# Define fixity of soil base nodes
for {set i 1} {$i<=[expr $ColWidthMesh+1]} {incr i 1} {
fix $i 0 1
}
# Define equal DOF for dashpot-soil nodes
equalDOF 1 4001 1
# Define equal DOF for base soil nodes
for {set i 2} {$i<=[expr $ColWidthMesh+1]} {incr i 1} {
equalDOF 1 $i 1
}
# Define equal DOF for boundary soil nodes
for {set i [expr $ColWidthMesh+2]} {$i<=[expr $NODESNUM-$ColWidthMesh]} {incr i [expr $ColWidthMesh+1]} {
for {set ind 2} {$ind<=$ColWidthMesh} {incr ind 1} {
equalDOF $i [expr $i+$ind] 1 2
}
}
puts "Finished creating all boundary conditions and equalDOF..."
#-------------------------------------------------------------------------------------------
# 5. DEFINE MATERIALS FOR SOIL ELEMENTS
#-------------------------------------------------------------------------------------------
for {set k 1} {$k <= $ColHeightMesh} {incr k 1} {
set gcurve "$rho $G($k) $bulk($k) $cohesion 0.1 $phi $refPress($k) $pressCoeff 20"; #AUTOMATIC GENERATION
nDMaterial PressureIndependMultiYield $k 2 {*}$gcurve
}
#-------------------------------------------------------------------------------------------
# 6. DEFINE SOIL ELEMENTS
#-------------------------------------------------------------------------------------------
set wgtX 0.0
set wgtY [expr -$g*$SOILDENSITY]
set EleTag 1
for {set j 0} {$j<=[expr $ColHeightMesh-1]} {incr j 1} {
for {set i 1} {$i <= $ColWidthMesh} {incr i 1} {
set index [expr $i+($ColWidthMesh+1)*$j]
set iNode $index
set jNode [expr $index+1]
set kNode [expr $index+$ColWidthMesh+2]
set lNode [expr $index+$ColWidthMesh+1]
element quad $EleTag $iNode $jNode $kNode $lNode 1.0 "PlaneStrain" [expr $j+1] 0.0 0.0 $wgtX $wgtY
set EleTag [expr $EleTag+1]
}
}
puts "Finished creating all soil elements..."
#-------------------------------------------------------------------------------------------
# 7. DEFINE MATERIAL AND ELEMENTS FOR VISCOUS DAMPERS
#-------------------------------------------------------------------------------------------
# Dashpot coefficient
set sizeEleX $ColWidth
set mC [expr $sizeEleX*$ROCKDENSITY*$ROCKVS]
# material
uniaxialMaterial Viscous 500 $mC 1
# elements
element zeroLength 5000 4000 4001 -mat 500 -dir 1
puts "Finished creating dashpot material and element..."
#-----------------------------------------------------------------------------------------------------------
# 9. APPLY GRAVITY LOADING
#-----------------------------------------------------------------------------------------------------------
constraints Transformation
test NormDispIncr 1e-5 100 0
algorithm Newton
numberer RCM
system ProfileSPD
integrator Newmark 0.5 0.25
analysis Transient
analyze 10 5.0e2
puts "Finished with elastic gravity analysis..."
# update materials to consider plastic behavior
for {set k 1} {$k <= $numLayers} {incr k} {
updateMaterialStage -material $k -stage 1
}
# plastic gravity loading
analyze 40 5.0e2
puts "Finished with plastic gravity analysis..."
#-------------------------------------------------------------------------------------------
# 8. DEFINE TIME SERIES
#-------------------------------------------------------------------------------------------
#---GROUND MOTION PARAMETERS
# time step in ground motion record
set motionDT $RECDT
# number of steps in ground motion record
set motionSteps $RECSTEPS
#---WAVELENGTH PARAMETERS
# damping coefficients (RAYLEIGH)
#set damp $SOILDAMPING
set damp 0.01
set a0 [expr 2*$damp*$omega1*$omega2/($omega1 + $omega2)]
set a1 [expr 2*$damp/($omega1 + $omega2)]
#---ANALYSIS PARAMETERS
# Newmark parameters
set gamma 0.5
set beta 0.25
#-------------------------------------------------------------------------------------------
# 9. CREATE POST-GRAVITY RECORDERS
#-------------------------------------------------------------------------------------------
# reset time and analysis
setTime 0.0
wipeAnalysis
# NOTE: RECORDERS SET FOR 2.00m column width
for {set k 0} {$k <= [expr $numLayers-1]} {incr k 1} {
set FFMaccname "OUTPUT_1D\\OS_temp_1D\\FFMacc$k.out"
set FFMdispname "OUTPUT_1D\\OS_temp_1D\\FFMdisp$k.out"
set ELEstrainname "OUTPUT_1D\\OS_temp_1D\\ELEstrain$k.out"
set ELEstressname "OUTPUT_1D\\OS_temp_1D\\ELEstress$k.out"
set Backbonename "OUTPUT_1D\\OS_temp_1D\\Backbone$k.out"
set FFMaccnode [expr 3*($numLayers+1)-(3*$k+1)]
set FFMdispnode $FFMaccnode
set ELEstrainele [expr 2*$numLayers-(2*$k+1)]
set ELEstressele $ELEstrainele
recorder Node -file $FFMaccname -time -dT $motionDT -node $FFMaccnode -dof 1 accel
recorder Node -file $FFMdispname -time -dT $motionDT -node $FFMdispnode -dof 1 accel
recorder Element -file $ELEstrainname -time -ele $ELEstrainele material 1 strain
recorder Element -file $ELEstressname -time -ele $ELEstressele material 1 stress
recorder Element -file $Backbonename -ele $ELEstressele material 1 backbone 80.0
}
puts "Finished creating all recorders..."
#-------------------------------------------------------------------------------------------
# 10. DYNAMIC ANALYSIS
#-------------------------------------------------------------------------------------------
# define constant factor for applied velocity
set cFactor [expr $sizeEleX*$ROCKDENSITY*$ROCKVS]
# define velocity time history file
set velocityFile Velocity1D.out
# timeseries object for applied force history
set timeSeries "Path -dt $motionDT -filePath $velocityFile -factor $cFactor "
# loading object
pattern Plain 1 $timeSeries {
load 1 1.0 0.0
}
puts "Dynamic loading created..."
#---DETERMINE STABLE ANALYSIS TIME STEP USING CFL CONDITION
# duration of ground motion (s)
set duration [expr $motionDT*$motionSteps]
# trial analysis time step
set kTrial [expr $sizeEleX/(pow($SOILVSff,0.5))]
# define time step and number of steps for analysis
if { $motionDT <= $kTrial } {
set nSteps $motionSteps
set dT $motionDT
} else {
set nSteps [expr floor($duration/$kTrial)+1]
set dT [expr $duration/$nSteps]
}
puts "number of steps in analysis: $nSteps"
puts "analysis time step: $dT"
# analysis objects
constraints Transformation
test NormDispIncr 1e-4 1000 0
algorithm Newton
numberer RCM
system ProfileSPD
integrator Newmark $gamma $beta
rayleigh $a0 $a1 0.0 0.0
analysis Transient
set ok [analyze $nSteps $dT]; # actually perform analysis; returns ok=0 if analysis was successful
puts "Finished with dynamic analysis..."
wipe all
quit
wipe all;
source INPUT_1D\\OPENSEES1DINPUT.txt
file mkdir OUTPUT_1D\\OS_temp_1D
#-------------------------------------------------------------------------------------------
# INPUT VARIABLES - PRELIMINARY CALCULATIONS
#-------------------------------------------------------------------------------------------
## Parameters that are sourced
#SOILHEIGHT
#SOILVS
#SOILDENSITY
#SOILPOISSON
#cohesion
#phi
#ROCKVS
#ROCKDENSITY
#RECDT
#RECSTEPS
#Define constants - calculate soil parameters - depth for output
set g 9.81
set pi 3.141592654
set thick 1.0; #Element thickness
set ColWidth 2
set WidthMesh 1
set ColWidthMesh [expr $ColWidth/$WidthMesh]
set HeightMesh 1
set ColHeightMesh [expr $SOILHEIGHT/$HeightMesh]
set pressCoeff 0.0
#########################################################
# #
# DEFINE LAYER PROPERTIES #
# #
#########################################################
set d 0.2
set numLayers $ColHeightMesh
set nu $SOILPOISSON
set rho $SOILDENSITY
set cohesion [expr $cohesion*1.0]
set SumVs 0.0
#Set Vsref to match Vs30=SOILVS
for {set k 1} {$k <= 30} {incr k 1} {
set SumVs [expr $SumVs+$SOILVS*pow(($k),$d/2)]
}
set Vs30 [expr $SumVs/30]
set VsRat [expr $SOILVS/$Vs30]
set SOILVSff [expr $SOILVS*$VsRat]
set profileFile [open OUTPUT_1D\\VsProfile.txt w]
# Vs for each layer (m/s)
for {set k $numLayers} {$k >= 1} {incr k -1} {
set Vs($k) [expr $SOILVSff*pow(($numLayers+1-$k),$d/2)]
puts $profileFile "$Vs($k)"
}
close $profileFile
# soil shear modulus for each layer (kPa)
for {set k 1} {$k <= $numLayers} {incr k 1} {
set G($k) [expr $rho*$Vs($k)*$Vs($k)]
}
# soil elastic modulus for each layer (kPa)
for {set k 1} {$k <= $numLayers} {incr k 1} {
set E($k) [expr 2*$G($k)*(1+$nu)]
}
# soil bulk modulus for each layer (kPa)
for {set k 1} {$k <= $numLayers} {incr k 1} {
set bulk($k) [expr $E($k)/(3*(1-2*$nu))]
}
# reference pressure for each layer (kPa)
for {set k 1} {$k <= $numLayers} {incr k 1} {
set refPress($k) 80.0; #Constant
}
# calculate soil period and frequencies
set SOILPERIOD 0.0
for {set k $numLayers} {$k >= 1} {incr k -1} {
set SOILPERIOD [expr $SOILPERIOD + $HeightMesh/$Vs($k)]
}
set SOILPERIOD [expr 4*$SOILPERIOD]
set omega1 [expr 2*$pi / $SOILPERIOD]
set omega2 [expr 5*$omega1]
#-------------------------------------------------------------------------------------------
# 1. DEFINE MODEL BUILDER
#-------------------------------------------------------------------------------------------
model basic -ndm 2 -ndf 2; # 2 spacial dimensions,2 DOF's per node
#----------------------------------------------------------------------------------------------------
# 2. DEFINE SOIL NODES
#----------------------------------------------------------------------------------------------------
set ntag 0
for {set iy 0} {$iy<=$ColHeightMesh} {incr iy 1} {
set nodey [expr $iy*$HeightMesh]
for {set ix 0} {$ix <= $ColWidthMesh} {incr ix 1} {
set nodex [expr $ix*$WidthMesh]
set ntag [expr $ntag+1]
node $ntag $nodex $nodey
}
}
puts "Finished creating soil nodes..."
set NODESNUM $ntag
#-----------------------------------------------------------------------------------------------------------
# 3. DEFINE DASHPOT NODES
#-----------------------------------------------------------------------------------------------------------
node 4000 0.0 0.0
node 4001 0.0 0.0
#-----------------------------------------------------------------------------------------------------------
# 4. DEFINE FIXITY AND EQUAL DOF
#-----------------------------------------------------------------------------------------------------------
# Define fixity of dashpot nodes
fix 4000 1 1
fix 4001 0 1
# Define fixity of soil base nodes
for {set i 1} {$i<=[expr $ColWidthMesh+1]} {incr i 1} {
fix $i 0 1
}
# Define equal DOF for dashpot-soil nodes
equalDOF 1 4001 1
# Define equal DOF for base soil nodes
for {set i 2} {$i<=[expr $ColWidthMesh+1]} {incr i 1} {
equalDOF 1 $i 1
}
# Define equal DOF for boundary soil nodes
for {set i [expr $ColWidthMesh+2]} {$i<=[expr $NODESNUM-$ColWidthMesh]} {incr i [expr $ColWidthMesh+1]} {
for {set ind 2} {$ind<=$ColWidthMesh} {incr ind 1} {
equalDOF $i [expr $i+$ind] 1 2
}
}
puts "Finished creating all boundary conditions and equalDOF..."
#-------------------------------------------------------------------------------------------
# 5. DEFINE MATERIALS FOR SOIL ELEMENTS
#-------------------------------------------------------------------------------------------
for {set k 1} {$k <= $ColHeightMesh} {incr k 1} {
set gcurve "$rho $G($k) $bulk($k) $cohesion 0.1 $phi $refPress($k) $pressCoeff 20"; #AUTOMATIC GENERATION
nDMaterial PressureIndependMultiYield $k 2 {*}$gcurve
}
#-------------------------------------------------------------------------------------------
# 6. DEFINE SOIL ELEMENTS
#-------------------------------------------------------------------------------------------
set wgtX 0.0
set wgtY [expr -$g*$SOILDENSITY]
set EleTag 1
for {set j 0} {$j<=[expr $ColHeightMesh-1]} {incr j 1} {
for {set i 1} {$i <= $ColWidthMesh} {incr i 1} {
set index [expr $i+($ColWidthMesh+1)*$j]
set iNode $index
set jNode [expr $index+1]
set kNode [expr $index+$ColWidthMesh+2]
set lNode [expr $index+$ColWidthMesh+1]
element quad $EleTag $iNode $jNode $kNode $lNode 1.0 "PlaneStrain" [expr $j+1] 0.0 0.0 $wgtX $wgtY
set EleTag [expr $EleTag+1]
}
}
puts "Finished creating all soil elements..."
#-------------------------------------------------------------------------------------------
# 7. DEFINE MATERIAL AND ELEMENTS FOR VISCOUS DAMPERS
#-------------------------------------------------------------------------------------------
# Dashpot coefficient
set sizeEleX $ColWidth
set mC [expr $sizeEleX*$ROCKDENSITY*$ROCKVS]
# material
uniaxialMaterial Viscous 500 $mC 1
# elements
element zeroLength 5000 4000 4001 -mat 500 -dir 1
puts "Finished creating dashpot material and element..."
#-----------------------------------------------------------------------------------------------------------
# 9. APPLY GRAVITY LOADING
#-----------------------------------------------------------------------------------------------------------
constraints Transformation
test NormDispIncr 1e-5 100 0
algorithm Newton
numberer RCM
system ProfileSPD
integrator Newmark 0.5 0.25
analysis Transient
analyze 10 5.0e2
puts "Finished with elastic gravity analysis..."
# update materials to consider plastic behavior
for {set k 1} {$k <= $numLayers} {incr k} {
updateMaterialStage -material $k -stage 1
}
# plastic gravity loading
analyze 40 5.0e2
puts "Finished with plastic gravity analysis..."
#-------------------------------------------------------------------------------------------
# 8. DEFINE TIME SERIES
#-------------------------------------------------------------------------------------------
#---GROUND MOTION PARAMETERS
# time step in ground motion record
set motionDT $RECDT
# number of steps in ground motion record
set motionSteps $RECSTEPS
#---WAVELENGTH PARAMETERS
# damping coefficients (RAYLEIGH)
#set damp $SOILDAMPING
set damp 0.01
set a0 [expr 2*$damp*$omega1*$omega2/($omega1 + $omega2)]
set a1 [expr 2*$damp/($omega1 + $omega2)]
#---ANALYSIS PARAMETERS
# Newmark parameters
set gamma 0.5
set beta 0.25
#-------------------------------------------------------------------------------------------
# 9. CREATE POST-GRAVITY RECORDERS
#-------------------------------------------------------------------------------------------
# reset time and analysis
setTime 0.0
wipeAnalysis
# NOTE: RECORDERS SET FOR 2.00m column width
for {set k 0} {$k <= [expr $numLayers-1]} {incr k 1} {
set FFMaccname "OUTPUT_1D\\OS_temp_1D\\FFMacc$k.out"
set FFMdispname "OUTPUT_1D\\OS_temp_1D\\FFMdisp$k.out"
set ELEstrainname "OUTPUT_1D\\OS_temp_1D\\ELEstrain$k.out"
set ELEstressname "OUTPUT_1D\\OS_temp_1D\\ELEstress$k.out"
set Backbonename "OUTPUT_1D\\OS_temp_1D\\Backbone$k.out"
set FFMaccnode [expr 3*($numLayers+1)-(3*$k+1)]
set FFMdispnode $FFMaccnode
set ELEstrainele [expr 2*$numLayers-(2*$k+1)]
set ELEstressele $ELEstrainele
recorder Node -file $FFMaccname -time -dT $motionDT -node $FFMaccnode -dof 1 accel
recorder Node -file $FFMdispname -time -dT $motionDT -node $FFMdispnode -dof 1 accel
recorder Element -file $ELEstrainname -time -ele $ELEstrainele material 1 strain
recorder Element -file $ELEstressname -time -ele $ELEstressele material 1 stress
recorder Element -file $Backbonename -ele $ELEstressele material 1 backbone 80.0
}
puts "Finished creating all recorders..."
#-------------------------------------------------------------------------------------------
# 10. DYNAMIC ANALYSIS
#-------------------------------------------------------------------------------------------
# define constant factor for applied velocity
set cFactor [expr $sizeEleX*$ROCKDENSITY*$ROCKVS]
# define velocity time history file
set velocityFile Velocity1D.out
# timeseries object for applied force history
set timeSeries "Path -dt $motionDT -filePath $velocityFile -factor $cFactor "
# loading object
pattern Plain 1 $timeSeries {
load 1 1.0 0.0
}
puts "Dynamic loading created..."
#---DETERMINE STABLE ANALYSIS TIME STEP USING CFL CONDITION
# duration of ground motion (s)
set duration [expr $motionDT*$motionSteps]
# trial analysis time step
set kTrial [expr $sizeEleX/(pow($SOILVSff,0.5))]
# define time step and number of steps for analysis
if { $motionDT <= $kTrial } {
set nSteps $motionSteps
set dT $motionDT
} else {
set nSteps [expr floor($duration/$kTrial)+1]
set dT [expr $duration/$nSteps]
}
puts "number of steps in analysis: $nSteps"
puts "analysis time step: $dT"
# analysis objects
constraints Transformation
test NormDispIncr 1e-4 1000 0
algorithm Newton
numberer RCM
system ProfileSPD
integrator Newmark $gamma $beta
rayleigh $a0 $a1 0.0 0.0
analysis Transient
set ok [analyze $nSteps $dT]; # actually perform analysis; returns ok=0 if analysis was successful
puts "Finished with dynamic analysis..."
wipe all
quit