I'm trying to construct a model of one block resting on top of another that I could then apply a lateral load to. I'm using quadrilateral elements and zeroLengthContact elements to connect the blocks. If I apply the lateral load in one load step it seems to work fine. However, if I attempt to apply a lateral load in more than one step, the first step works fine but any step beyond the first produces an instability. The displacements essentially go to infinity, as if OpenSees does not recognize that the two blocks are still contacting after the first load step.
Any thoughts as to why this may be occurring? My goal is to ultimately perform a dynamic analysis of these blocks therefore I'm going to need to find a better way to load them than applying single step lateral loads and resetting the model after each step.
#####################
Here is the code to construct the model:
######################
wipe
model BasicBuilder -ndm 2 -ndf 2
# Define Nodes
set H 12; #inches
set L 12; #inches
set D 12; #inches
node 1 0 0
node 2 $L 0
node 3 $L $H
node 4 0 $H
node 5 0 $H
node 6 $L $H
node 7 $L [expr 2*$H]
node 8 0 [expr 2*$H]
# Define Fixities
fix 1 1 1
fix 2 1 1
nDMaterial ElasticIsotropic 1 6000 0.25; # E (ksi)
# Define Elements
set rho 0.000939236; # density, kips/cubic inch
set bulkmod 3000; # Bulk modulus, ksi
set V [expr $H*$L*$D]; # Volume of element
set k 0.1; # penalty factor
set T [expr $k*$bulkmod*$L*$D*$L*$D/$V]
set P [expr ($rho*$L*$H*$D)/2]
set g 32.2
set m [expr 1000*$rho*$V/$g]
element quad 1 1 2 3 4 $D PlaneStrain 1
element quad 2 5 6 7 8 $D PlaneStrain 1
element zeroLengthContact2D 3 5 4 $T $T 0.6 -normal 0 1
element zeroLengthContact2D 4 6 3 $T $T 0.6 -normal 0 1
mass 5 $m 0
mass 6 $m 0
# Set Axial Load
set axialLoad -$P
# Define Constant Axial Load
pattern Plain 1 Constant {
load 5 0 $axialLoad
load 6 0 $axialLoad
};
# Define Analysis Paramters
system BandGeneral
integrator LoadControl 1
test NormUnbalance 1.0e-6 100
numberer Plain
constraints Plain
algorithm Linear
analysis Static
# Perform 1 analysis for constant axial load
analyze 1
loadConst -time 0.0
#######################
Here is the code to apply the lateral load:
#######################
wipe
source BuildModelSimple.tcl
# Define Lateral Load Pattern
set P 100
pattern Plain 2 Linear {
load 5 $P 0
load 6 $P 0
}
# Recorders
recorder Node -file NodeDisplacementsSimple.dat -node 1 2 3 4 5 6 7 8 -dof 1 2 disp
# Perform Push-Over
constraints Plain
numberer Plain
system BandGeneral
test EnergyIncr 1.0e-9 100
algorithm Linear
integrator LoadControl 1
analysis Static
analyze 1
Note: LoadControl 1 and analyze 1 works fine, but if I change LoadControl to 0.1 and analyze to 10, the displacements diverge on the second load step.
Instabilities with zeroLengthContact
Moderators: silvia, selimgunay, Moderators
The example pushes two blocks with applied axial load = 0.8, and lateral load=100, with frictional coefficient=0.6 between blocks.
The static solution is non-existent since the system cannot stay in equilibrium.
My suggestion:
(1) reduce lateral load to a reasonable number
(2) set k as a large penalty factor. The penalty is used to reinforce the frictional constraint, it has to be a reasonably large number. The larger the better, but too large will cause numerical instability.
(3) change to transient solver to solve a dynamic problem. The block will accelerate.
The static solution is non-existent since the system cannot stay in equilibrium.
My suggestion:
(1) reduce lateral load to a reasonable number
(2) set k as a large penalty factor. The penalty is used to reinforce the frictional constraint, it has to be a reasonably large number. The larger the better, but too large will cause numerical instability.
(3) change to transient solver to solve a dynamic problem. The block will accelerate.
Hi Gang,
Thanks so much for your reply. My goal is to perform dynamic time history analyses of a structure modeled with only quadrilateral elements with zeroLengthContact elements connecting them. I have been able to get my 2 block system to work under a dynamic analysis, however even for a large scale earthquake I'm not seeing much displacement (roughly 0.005 inches). Could you take a look at my code and see if anything jumps out at you as being incorrect? The analysis is very sensitive to the size of the masses, if I double them I instantly get very large displacements (roughly 4 inches).
Thanks again for the help.
########################
Build the Model
########################
wipe
model BasicBuilder -ndm 2 -ndf 2
# Define Nodes
set H 12; #inches
set L 12; #inches
set D 12; #inches
node 1 0 0
node 2 $L 0
node 3 $L $H
node 4 0 $H
node 5 0 $H
node 6 $L $H
node 7 $L [expr 2*$H]
node 8 0 [expr 2*$H]
# Define Fixities
fix 1 1 1
fix 2 1 1
nDMaterial ElasticIsotropic 1 6E6 0.25; # E (psi)
# Define Elements
set rho 0.09375; # density, pounds/cubic inch
set bulkmod 3E6; # Bulk modulus, psi
set V [expr $H*$L*$D]; # Volume of element
set k 0.1; # penalty factor
set A [expr $H*$L]
set T [expr $k*$bulkmod*$A*$A/$V]
set P [expr ($rho*$V)/2]
set g 386.4
set m [expr $rho*$V/$g/2]
element quad 1 1 2 3 4 $D PlaneStrain 1
element quad 2 5 6 7 8 $D PlaneStrain 1
element zeroLengthContact2D 3 5 4 $T $T 0.6 -normal 0 1
element zeroLengthContact2D 4 6 3 $T $T 0.6 -normal 0 1
mass 5 $m 0
mass 6 $m 0
# Set Axial Load
set axialLoad -$P
# Define Constant Axial Load
pattern Plain 1 Constant {
load 5 0 $axialLoad
load 6 0 $axialLoad
};
# Define Analysis Paramters
system BandGeneral
integrator LoadControl 1
test NormUnbalance 1.0e-6 100
numberer Plain
constraints Plain
algorithm Linear
analysis Static
# Perform 1 analysis for constant axial load
analyze 1
loadConst -time 0.0
########################
Dynamic Analysis
########################
wipe
# Assemble Model and Apply Gravity
source BuildModelSimple.tcl
source ReadSMDFile.tcl
# Uniform Earthquake ground motion (uniform acceleration input at all support nodes)
set GMdirection 1; # ground-motion direction
set GMfile "TAK090" ; # ground-motion filenames
set GMfact 1; # ground-motion scaling factor
set g 386.4
# set up ground-motion-analysis parameters
set DtAnalysis [expr 0.01]; # time-step Dt for lateral analysis
set TmaxAnalysis [expr 40]; # maximum duration of ground-motion analysis -- should be 50*$sec
# Define Analysis Paramters
constraints Penalty 1E20 1E20
numberer RCM
system UmfPack
test EnergyIncr 1e-08 1000
algorithm NewtonLineSearch 0.6
integrator Newmark 0.5 0.25
analysis Transient
# RAYLEIGH damping parameters, Where to put M/K-prop damping, switches (http://opensees.berkeley.edu/OpenSees/m ... l/1099.htm)
# D=$alphaM*M + $betaKcurr*Kcurrent + $betaKcomm*KlastCommit + $beatKinit*$Kinitial
set xDamp 0.02; # damping ratio
set MpropSwitch 1.0;
set KcurrSwitch 0.0;
set KcommSwitch 0.0;
set KinitSwitch 1.0;
set nEigenI 1; # mode 1
set nEigenJ 2; # mode 2
set lambdaN [eigen [expr $nEigenJ]]; # eigenvalue analysis for nEigenJ modes
set lambdaI [lindex $lambdaN [expr $nEigenI-1]]; # eigenvalue mode i
set lambdaJ [lindex $lambdaN [expr $nEigenJ-1]]; # eigenvalue mode j
set omegaI [expr pow($lambdaI,0.5)];
set omegaJ [expr pow($lambdaJ,0.5)];
set alphaM [expr $MpropSwitch*$xDamp*(2*$omegaI*$omegaJ)/($omegaI+$omegaJ)]; # M-prop. damping; D = alphaM*M
set betaKcurr [expr $KcurrSwitch*2.*$xDamp/($omegaI+$omegaJ)]; # current-K; +beatKcurr*KCurrent
set betaKcomm [expr $KcommSwitch*2.*$xDamp/($omegaI+$omegaJ)]; # last-committed K; +betaKcomm*KlastCommitt
set betaKinit [expr $KinitSwitch*2.*$xDamp/($omegaI+$omegaJ)]; # initial-K; +beatKinit*Kini
rayleigh $alphaM $betaKcurr $betaKinit $betaKcomm; # RAYLEIGH damping
# --------------------------------- perform Dynamic Ground-Motion Analysis
# the following commands are unique to the Uniform Earthquake excitation
set IDloadTag 400; # for uniformSupport excitation
# Uniform EXCITATION: acceleration input
set inFile $GMfile.at2
set outFile $GMfile.g3; # set variable holding new filename (PEER files have .at2/dt2 extension)
ReadSMDFile $inFile $outFile dt; # call procedure to convert the ground-motion file
set GMfatt [expr $g*$GMfact]; # data in input file is in g Unifts -- ACCELERATION TH
set AccelSeries "Series -dt $dt -filePath $outFile -factor $GMfatt"; # time series information
pattern UniformExcitation $IDloadTag $GMdirection -accel $AccelSeries ; # create Unifform excitation
# Recorders
recorder Node -file DynamicDisplacements.tcl -node 5 6 -dof 1 2 disp
recorder Node -file DynamicReactions.tcl -node 1 2 -dof 2 reaction
set Nsteps [expr int($TmaxAnalysis/$DtAnalysis)];
set ok [analyze $Nsteps $DtAnalysis]; # actually perform analysis; returns ok=0 if analysis was successful
puts "Ground Motion Done. End Time: [getTime]"
Thanks so much for your reply. My goal is to perform dynamic time history analyses of a structure modeled with only quadrilateral elements with zeroLengthContact elements connecting them. I have been able to get my 2 block system to work under a dynamic analysis, however even for a large scale earthquake I'm not seeing much displacement (roughly 0.005 inches). Could you take a look at my code and see if anything jumps out at you as being incorrect? The analysis is very sensitive to the size of the masses, if I double them I instantly get very large displacements (roughly 4 inches).
Thanks again for the help.
########################
Build the Model
########################
wipe
model BasicBuilder -ndm 2 -ndf 2
# Define Nodes
set H 12; #inches
set L 12; #inches
set D 12; #inches
node 1 0 0
node 2 $L 0
node 3 $L $H
node 4 0 $H
node 5 0 $H
node 6 $L $H
node 7 $L [expr 2*$H]
node 8 0 [expr 2*$H]
# Define Fixities
fix 1 1 1
fix 2 1 1
nDMaterial ElasticIsotropic 1 6E6 0.25; # E (psi)
# Define Elements
set rho 0.09375; # density, pounds/cubic inch
set bulkmod 3E6; # Bulk modulus, psi
set V [expr $H*$L*$D]; # Volume of element
set k 0.1; # penalty factor
set A [expr $H*$L]
set T [expr $k*$bulkmod*$A*$A/$V]
set P [expr ($rho*$V)/2]
set g 386.4
set m [expr $rho*$V/$g/2]
element quad 1 1 2 3 4 $D PlaneStrain 1
element quad 2 5 6 7 8 $D PlaneStrain 1
element zeroLengthContact2D 3 5 4 $T $T 0.6 -normal 0 1
element zeroLengthContact2D 4 6 3 $T $T 0.6 -normal 0 1
mass 5 $m 0
mass 6 $m 0
# Set Axial Load
set axialLoad -$P
# Define Constant Axial Load
pattern Plain 1 Constant {
load 5 0 $axialLoad
load 6 0 $axialLoad
};
# Define Analysis Paramters
system BandGeneral
integrator LoadControl 1
test NormUnbalance 1.0e-6 100
numberer Plain
constraints Plain
algorithm Linear
analysis Static
# Perform 1 analysis for constant axial load
analyze 1
loadConst -time 0.0
########################
Dynamic Analysis
########################
wipe
# Assemble Model and Apply Gravity
source BuildModelSimple.tcl
source ReadSMDFile.tcl
# Uniform Earthquake ground motion (uniform acceleration input at all support nodes)
set GMdirection 1; # ground-motion direction
set GMfile "TAK090" ; # ground-motion filenames
set GMfact 1; # ground-motion scaling factor
set g 386.4
# set up ground-motion-analysis parameters
set DtAnalysis [expr 0.01]; # time-step Dt for lateral analysis
set TmaxAnalysis [expr 40]; # maximum duration of ground-motion analysis -- should be 50*$sec
# Define Analysis Paramters
constraints Penalty 1E20 1E20
numberer RCM
system UmfPack
test EnergyIncr 1e-08 1000
algorithm NewtonLineSearch 0.6
integrator Newmark 0.5 0.25
analysis Transient
# RAYLEIGH damping parameters, Where to put M/K-prop damping, switches (http://opensees.berkeley.edu/OpenSees/m ... l/1099.htm)
# D=$alphaM*M + $betaKcurr*Kcurrent + $betaKcomm*KlastCommit + $beatKinit*$Kinitial
set xDamp 0.02; # damping ratio
set MpropSwitch 1.0;
set KcurrSwitch 0.0;
set KcommSwitch 0.0;
set KinitSwitch 1.0;
set nEigenI 1; # mode 1
set nEigenJ 2; # mode 2
set lambdaN [eigen [expr $nEigenJ]]; # eigenvalue analysis for nEigenJ modes
set lambdaI [lindex $lambdaN [expr $nEigenI-1]]; # eigenvalue mode i
set lambdaJ [lindex $lambdaN [expr $nEigenJ-1]]; # eigenvalue mode j
set omegaI [expr pow($lambdaI,0.5)];
set omegaJ [expr pow($lambdaJ,0.5)];
set alphaM [expr $MpropSwitch*$xDamp*(2*$omegaI*$omegaJ)/($omegaI+$omegaJ)]; # M-prop. damping; D = alphaM*M
set betaKcurr [expr $KcurrSwitch*2.*$xDamp/($omegaI+$omegaJ)]; # current-K; +beatKcurr*KCurrent
set betaKcomm [expr $KcommSwitch*2.*$xDamp/($omegaI+$omegaJ)]; # last-committed K; +betaKcomm*KlastCommitt
set betaKinit [expr $KinitSwitch*2.*$xDamp/($omegaI+$omegaJ)]; # initial-K; +beatKinit*Kini
rayleigh $alphaM $betaKcurr $betaKinit $betaKcomm; # RAYLEIGH damping
# --------------------------------- perform Dynamic Ground-Motion Analysis
# the following commands are unique to the Uniform Earthquake excitation
set IDloadTag 400; # for uniformSupport excitation
# Uniform EXCITATION: acceleration input
set inFile $GMfile.at2
set outFile $GMfile.g3; # set variable holding new filename (PEER files have .at2/dt2 extension)
ReadSMDFile $inFile $outFile dt; # call procedure to convert the ground-motion file
set GMfatt [expr $g*$GMfact]; # data in input file is in g Unifts -- ACCELERATION TH
set AccelSeries "Series -dt $dt -filePath $outFile -factor $GMfatt"; # time series information
pattern UniformExcitation $IDloadTag $GMdirection -accel $AccelSeries ; # create Unifform excitation
# Recorders
recorder Node -file DynamicDisplacements.tcl -node 5 6 -dof 1 2 disp
recorder Node -file DynamicReactions.tcl -node 1 2 -dof 2 reaction
set Nsteps [expr int($TmaxAnalysis/$DtAnalysis)];
set ok [analyze $Nsteps $DtAnalysis]; # actually perform analysis; returns ok=0 if analysis was successful
puts "Ground Motion Done. End Time: [getTime]"
You need to set penalty a much larger number to reinforce contact law.
for example, please change the following line:
set k 10000; # penalty factor
If you record the displacement, you can check by hand at what time, and at what acceleration that the sliding starts, theoretically.
I cannot test your code since they are not complete. Please include everything, including BuildModelSimple.tcl and
ReadSMDFile.tcl and ground motion to send to me via gwang@ust.hk, so I can test for you.
for example, please change the following line:
set k 10000; # penalty factor
If you record the displacement, you can check by hand at what time, and at what acceleration that the sliding starts, theoretically.
I cannot test your code since they are not complete. Please include everything, including BuildModelSimple.tcl and
ReadSMDFile.tcl and ground motion to send to me via gwang@ust.hk, so I can test for you.