Hi all,
I have created a model of a pile (elasitcBeamColumn elements) in liquefiable sands. The model is resting on a be of L-K dashpots, and the pile and soil are connected by small elasticBeamColumn elements. The ground excitation applied is a cyclic sinusoidal motion.
The model runs the first 3 analyses (gravity loads for soil elastic, plastic, and gravity loads with pile added). And it will partially run the fourth analysis (ground excitation). The earthquake is 20sec long with time steps of 0.01sec - 2000 time steps in total. The model will easily converge to a solution for the first 20-40 steps (depending on integrator, algorithm, constraints etc that are being used) then will actually diverge, eventually crashing OpenSees.
Has anyone come across this situation before (and found the solution)? I have checked the output for the steps that do converge and nothing is flying off into space - all the responses and deformations are as they should be. Neither does the onset of divergence correspond with the change in direction of earthquake shaking.
Could someone please help me, I'm at a loss.
The code follows in case anyone is interested:
#################################################
#
# Program models a pushover case in liquefiable sand employing soil-pile interface slip and gapping (modelled by a
# weaker zone of elements on the outside of the pile
# Gobal co-ordinate system: y +ve upwards, x +ve to the right, z +ve into the page
# Groundwater table at 3m depth
#
#################################################
# Unit system in this file: SI
set g 9.81
wipe
##########
# Set up soil elements
##########
model basic -ndm 3 -ndf 3
# add Lysmer-Kuhlemeyer damping for soil
set shearcoeff 564
set normcoeff 564
uniaxialMaterial Viscous 4 $normcoeff 1
uniaxialMaterial Viscous 5 $shearcoeff 1
set pi [expr 4*atan(1.0)]
set deg [expr $pi/180.]
# create normal direction dashpots
node 325 0.0 -9.0 0.5
node 396 -0.353 -9.0 0.353
node 326 -0.5 -9.0 0.0
node 327 -2.5 -9.0 0.0
node 328 -2.5 -9.0 2.5
node 329 0.0 -9.0 2.5
node 330 2.5 -9.0 2.5
node 331 2.5 -9.0 0.0
node 332 0.5 -9.0 0.0
node 397 0.353 -9.0 0.353
node 366 -0.29 -9.0 0.0
node 3107 [expr -0.29*sin(45*$deg)] -9.0 [expr 0.29*sin(45*$deg)]
node 367 0.0 -9.0 0.29
node 368 0.29 -9.0 0.0
node 3106 [expr 0.29*sin(45*$deg)] -9.0 [expr 0.29*sin(45*$deg)]
node 625 0.0 -9.0 0.5
node 696 -0.353 -9.0 0.353
node 626 -0.5 -9.0 0.0
node 627 -2.5 -9.0 0.0
node 628 -2.5 -9.0 2.5
node 629 0.0 -9.0 2.5
node 630 2.5 -9.0 2.5
node 631 2.5 -9.0 0.0
node 632 0.5 -9.0 0.0
node 697 0.353 -9.0 0.353
node 666 -0.29 -9.0 0.0
node 6107 [expr -0.29*sin(45*$deg)] -9.0 [expr 0.29*sin(45*$deg)]
node 667 0.0 -9.0 0.29
node 668 0.29 -9.0 0.0
node 6106 [expr 0.29*sin(45*$deg)] -9.0 [expr 0.29*sin(45*$deg)]
# create shear direction dashpots
node 425 0.0 -9.0 0.5
node 496 -0.353 -9.0 0.353
node 426 -0.5 -9.0 0.0
node 427 -2.5 -9.0 0.0
node 428 -2.5 -9.0 2.5
node 429 0.0 -9.0 2.5
node 430 2.5 -9.0 2.5
node 431 2.5 -9.0 0.0
node 432 0.5 -9.0 0.0
node 497 0.353 -9.0 0.353
node 466 -0.29 -9.0 0.0
node 4107 [expr -0.29*sin(45*$deg)] -9.0 [expr 0.29*sin(45*$deg)]
node 467 0.0 -9.0 0.29
node 468 0.29 -9.0 0.0
node 4106 [expr 0.29*sin(45*$deg)] -9.0 [expr 0.29*sin(45*$deg)]
node 725 0.0 -9.0 0.5
node 796 -0.353 -9.0 0.353
node 726 -0.5 -9.0 0.0
node 727 -2.5 -9.0 0.0
node 728 -2.5 -9.0 2.5
node 729 0.0 -9.0 2.5
node 730 2.5 -9.0 2.5
node 731 2.5 -9.0 0.0
node 732 0.5 -9.0 0.0
node 797 0.353 -9.0 0.353
node 766 -0.29 -9.0 0.0
node 7107 [expr -0.29*sin(45*$deg)] -9.0 [expr 0.29*sin(45*$deg)]
node 767 0.0 -9.0 0.29
node 768 0.29 -9.0 0.0
node 7106 [expr 0.29*sin(45*$deg)] -9.0 [expr 0.29*sin(45*$deg)]
fix 325 1 1 1
fix 396 1 1 1
fix 326 1 1 1
fix 327 1 1 1
fix 328 1 1 1
fix 329 1 1 1
fix 330 1 1 1
fix 331 1 1 1
fix 332 1 1 1
fix 397 1 1 1
fix 366 1 1 1
fix 3107 1 1 1
fix 367 1 1 1
fix 368 1 1 1
fix 3106 1 1 1
fix 625 1 0 1
fix 696 1 0 1
fix 626 1 0 1
fix 627 1 0 1
fix 628 1 0 1
fix 629 1 0 1
fix 630 1 0 1
fix 631 1 0 1
fix 632 1 0 1
fix 697 1 0 1
fix 666 1 0 1
fix 6107 1 0 1
fix 667 1 0 1
fix 668 1 0 1
fix 6106 1 0 1
fix 425 1 1 1
fix 496 1 1 1
fix 426 1 1 1
fix 427 1 1 1
fix 428 1 1 1
fix 429 1 1 1
fix 430 1 1 1
fix 431 1 1 1
fix 432 1 1 1
fix 497 1 1 1
fix 466 1 1 1
fix 4107 1 1 1
fix 467 1 1 1
fix 468 1 1 1
fix 4106 1 1 1
fix 725 0 1 1
fix 796 0 1 1
fix 726 0 1 1
fix 727 0 1 1
fix 728 0 1 1
fix 729 0 1 1
fix 730 0 1 1
fix 731 0 1 1
fix 732 0 1 1
fix 797 0 1 1
fix 766 0 1 1
fix 7107 0 1 1
fix 767 0 1 1
fix 768 0 1 1
fix 7106 0 1 1
# create dashpots
element zeroLength 300 325 625 -mat 4 -dir 2
element zeroLength 301 396 696 -mat 4 -dir 2
element zeroLength 302 326 626 -mat 4 -dir 2
element zeroLength 303 327 627 -mat 4 -dir 2
element zeroLength 304 328 628 -mat 4 -dir 2
element zeroLength 305 329 629 -mat 4 -dir 2
element zeroLength 306 330 630 -mat 4 -dir 2
element zeroLength 307 331 631 -mat 4 -dir 2
element zeroLength 308 332 632 -mat 4 -dir 2
element zeroLength 309 397 697 -mat 4 -dir 2
element zeroLength 310 366 666 -mat 4 -dir 2
element zeroLength 311 3107 6107 -mat 4 -dir 2
element zeroLength 312 367 667 -mat 4 -dir 2
element zeroLength 313 368 668 -mat 4 -dir 2
element zeroLength 314 3106 6106 -mat 4 -dir 2
element zeroLength 700 425 725 -mat 5 -dir 1
element zeroLength 701 496 796 -mat 5 -dir 1
element zeroLength 702 426 726 -mat 5 -dir 1
element zeroLength 703 427 727 -mat 5 -dir 1
element zeroLength 704 428 728 -mat 5 -dir 1
element zeroLength 705 429 729 -mat 5 -dir 1
element zeroLength 706 430 730 -mat 5 -dir 1
element zeroLength 707 431 731 -mat 5 -dir 1
element zeroLength 708 432 732 -mat 5 -dir 1
element zeroLength 709 497 797 -mat 5 -dir 1
element zeroLength 710 466 766 -mat 5 -dir 1
element zeroLength 711 4107 7107 -mat 5 -dir 1
element zeroLength 712 467 767 -mat 5 -dir 1
element zeroLength 713 468 768 -mat 5 -dir 1
element zeroLength 714 4106 7106 -mat 5 -dir 1
wipeAnalysis
model basic -ndm 3 -ndf 4
#define nodes
node 1 0.0 0.0 0.5
node 2 -0.5 0.0 0.0
node 3 -2.5 0.0 0.0
node 4 -2.5 0.0 2.5
node 5 0.0 0.0 2.5
node 6 2.5 0.0 2.5
node 7 2.5 0.0 0.0
node 8 0.5 0.0 0.0
node 9 0.0 -3.0 0.5
node 10 -0.5 -3.0 0.0
node 11 -2.5 -3.0 0.0
node 12 -2.5 -3.0 2.5
node 13 0.0 -3.0 2.5
node 14 2.5 -3.0 2.5
node 15 2.5 -3.0 0.0
node 16 0.5 -3.0 0.0
node 17 0.0 -6.0 0.5
node 18 -0.5 -6.0 0.0
node 19 -2.5 -6.0 0.0
node 20 -2.5 -6.0 2.5
node 21 0.0 -6.0 2.5
node 22 2.5 -6.0 2.5
node 23 2.5 -6.0 0.0
node 24 0.5 -6.0 0.0
node 25 0.0 -9.0 0.5
node 26 -0.5 -9.0 0.0
node 27 -2.5 -9.0 0.0
node 28 -2.5 -9.0 2.5
node 29 0.0 -9.0 2.5
node 30 2.5 -9.0 2.5
node 31 2.5 -9.0 0.0
node 32 0.5 -9.0 0.0
node 90 -0.353 0.0 0.353
node 91 0.353 0.0 0.353
node 92 -0.353 -3.0 0.353
node 93 0.353 -3.0 0.353
node 94 -0.353 -6.0 0.353
node 95 0.353 -6.0 0.353
node 96 -0.353 -9.0 0.353
node 97 0.353 -9.0 0.353
#contact material nodes
node 57 -0.29 0.0 0.0
node 58 0.0 0.0 0.29
node 59 0.29 0.0 0.0
node 60 -0.29 -3.0 0.0
node 61 0.0 -3.0 0.29
node 62 0.29 -3.0 0.0
node 63 -0.29 -6.0 0.0
node 64 0.0 -6.0 0.29
node 65 0.29 -6.0 0.0
node 66 -0.29 -9.0 0.0
node 67 0.0 -9.0 0.29
node 68 0.29 -9.0 0.0
# create nodes at soil-pile interface exactly 0.29m from centerline so nodes can have contact
set pi [expr 4*atan(1.0)]
set deg [expr $pi/180.]
node 100 [expr 0.29*sin(45*$deg)] 0.0 [expr 0.29*sin(45*$deg)]
node 101 [expr -0.29*sin(45*$deg)] 0.0 [expr 0.29*sin(45*$deg)]
node 102 [expr 0.29*sin(45*$deg)] -3.0 [expr 0.29*sin(45*$deg)]
node 103 [expr -0.29*sin(45*$deg)] -3.0 [expr 0.29*sin(45*$deg)]
node 104 [expr 0.29*sin(45*$deg)] -6.0 [expr 0.29*sin(45*$deg)]
node 105 [expr -0.29*sin(45*$deg)] -6.0 [expr 0.29*sin(45*$deg)]
node 106 [expr 0.29*sin(45*$deg)] -9.0 [expr 0.29*sin(45*$deg)]
node 107 [expr -0.29*sin(45*$deg)] -9.0 [expr 0.29*sin(45*$deg)]
# fix base of model
fix 25 0 0 1 0
fix 26 0 0 1 0
fix 27 0 0 1 0
fix 28 0 0 1 0
fix 29 0 0 1 0
fix 30 0 0 1 0
fix 31 0 0 1 0
fix 32 0 0 1 0
fix 96 0 0 1 0
fix 97 0 0 1 0
fix 66 0 0 1 0
fix 67 0 0 1 0
fix 68 0 0 1 0
fix 107 0 0 1 0
fix 106 0 0 1 0
# fix base of model to L-K dasphots
equalDOF 625 25 1 2 3
equalDOF 696 96 1 2 3
equalDOF 626 26 1 2 3
equalDOF 627 27 1 2 3
equalDOF 628 28 1 2 3
equalDOF 629 29 1 2 3
equalDOF 630 30 1 2 3
equalDOF 631 31 1 2 3
equalDOF 632 32 1 2 3
equalDOF 697 97 1 2 3
equalDOF 666 66 1 2 3
equalDOF 6107 107 1 2 3
equalDOF 667 67 1 2 3
equalDOF 668 68 1 2 3
equalDOF 6106 106 1 2 3
equalDOF 725 25 1 2 3
equalDOF 796 96 1 2 3
equalDOF 726 26 1 2 3
equalDOF 727 27 1 2 3
equalDOF 728 28 1 2 3
equalDOF 729 29 1 2 3
equalDOF 730 30 1 2 3
equalDOF 731 31 1 2 3
equalDOF 732 32 1 2 3
equalDOF 797 97 1 2 3
equalDOF 766 66 1 2 3
equalDOF 7107 107 1 2 3
equalDOF 767 67 1 2 3
equalDOF 768 68 1 2 3
equalDOF 7106 106 1 2 3
#equalDOF 625 25 2
#equalDOF 696 96 2
#equalDOF 626 26 2
#equalDOF 627 27 2
#equalDOF 628 28 2
#equalDOF 629 29 2
#equalDOF 630 30 2
#equalDOF 631 31 2
#equalDOF 632 32 2
#equalDOF 697 97 2
#equalDOF 666 66 2
#equalDOF 6107 107 2
#equalDOF 667 67 2
#equalDOF 668 68 2
#equalDOF 6106 106 2
#equalDOF 725 25 1
#equalDOF 796 96 1
#equalDOF 726 26 1
#equalDOF 727 27 1
#equalDOF 728 28 1
#equalDOF 729 29 1
#equalDOF 730 30 1
#equalDOF 731 31 1
#equalDOF 732 32 1
#equalDOF 797 97 1
#equalDOF 766 66 1
#equalDOF 7107 107 1
#equalDOF 767 67 1
#equalDOF 768 68 1
#equalDOF 7106 106 1
# fixes nodes on cross-section face from z (out of plane) movement on line of symmetry
fixZ [expr 0.0] 0 0 1 1 -tol 1e-4
fixZ [expr 2.5] 0 0 1 1 -tol 1e-4
# boundary conditions nodes slaved together in 1 (x disp) and 2 (y disp)
# creates periodic boundary conditions
equalDOF 3 7 1 2
equalDOF 11 15 1 2
equalDOF 19 23 1 2
equalDOF 27 31 1 2
equalDOF 4 6 1 2
equalDOF 12 14 1 2
equalDOF 20 22 1 2
equalDOF 28 30 1 2
# fix pore water condition at ground surface
fix 1 0 0 0 1
fix 2 0 0 0 1
fix 3 0 0 0 1
fix 4 0 0 0 1
fix 5 0 0 0 1
fix 6 0 0 0 1
fix 7 0 0 0 1
fix 8 0 0 0 1
fix 90 0 0 0 1
fix 91 0 0 0 1
fix 57 0 0 0 1
fix 101 0 0 0 1
fix 58 0 0 0 1
fix 100 0 0 0 1
fix 59 0 0 0 1
# fix pore water condition at 3m depth (depth of GWT)
fix 9 0 0 0 1
fix 10 0 0 0 1
fix 11 0 0 0 1
fix 12 0 0 0 1
fix 13 0 0 0 1
fix 14 0 0 0 1
fix 15 0 0 0 1
fix 16 0 0 0 1
fix 92 0 0 0 1
fix 93 0 0 0 1
fix 60 0 0 0 1
fix 103 0 0 0 1
fix 61 0 0 0 1
fix 102 0 0 0 1
fix 62 0 0 0 1
# soil material
nDMaterial PressureDependMultiYield 1 3 1.7 5.5e+004 1.5e+005 29 0.1 80 0.5 29 0.21 0 0 10 0.02 1
#nDMaterial PressureIndependMultiYield 1 3 1.3 1.3e+4 6.5e+04 18 0.1
# interface material
#nDMaterial PressureIndependMultiYield 9 3 1.7 7000 1.5e+005 0 0.1 17
#nDMaterial PressureIndependMultiYield 9 3 1.7 70000 1.5e+005 10 0.1 17
nDMaterial PressureDependMultiYield 9 3 1.7 7000 1.5e+005 17 0.1 80 0.5 17 0 0 0 0 0 0
set gravX 0
set gravZ 0
set gravY [expr -1*$g]
set xperm 1
set yperm 1
set zperm 1
set bulk 2.2000e+006
set fmass 1
# main body
element brickUP 1 3 4 90 2 11 12 92 10 1 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 2 11 12 92 10 19 20 94 18 1 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 3 19 20 94 18 27 28 96 26 1 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 4 90 4 5 1 92 12 13 9 1 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 5 92 12 13 9 94 20 21 17 1 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 6 94 20 21 17 96 28 29 25 1 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 7 1 5 6 91 9 13 14 93 1 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 8 9 13 14 93 17 21 22 95 1 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 9 17 21 22 95 25 29 30 97 1 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 10 91 6 7 8 93 14 15 16 1 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 11 93 14 15 16 95 22 23 24 1 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 12 95 22 23 24 97 30 31 32 1 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
# refined mesh around pile area
element brickUP 13 57 2 90 101 60 10 92 103 9 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 14 60 10 92 103 63 18 94 105 9 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 15 63 18 94 105 66 26 96 107 9 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 16 101 90 1 58 103 92 9 61 9 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 17 103 92 9 61 105 94 17 64 9 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 18 105 94 17 64 107 96 25 67 9 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 19 58 1 91 100 61 9 93 102 9 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 20 61 9 93 102 64 17 95 104 9 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 21 64 17 95 104 67 25 97 106 9 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 22 100 91 8 59 102 93 16 62 9 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 23 102 93 16 62 104 95 24 65 9 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
element brickUP 24 104 95 24 65 106 97 32 68 9 $bulk $fmass $xperm $yperm $zperm $gravX $gravY $gravZ
set NumElements 24
set NodeList {}
for {set i 1} {$i<=32} {incr i 1} {
lappend NodeList [expr $i]
}
for {set i 90} {$i<=97} {incr i 1} {
lappend NodeList [expr $i]
}
for {set i 57} {$i<=68} {incr i 1} {
lappend NodeList [expr $i]
}
for {set i 100} {$i<=107} {incr i 1} {
lappend NodeList [expr $i]
}
set SoilEleList {1 2 3 4 5 6 7 8 9 10 11 12}
set contactEleList {13 14 15 16 17 18 19 20 21 22 23 24}
########################
# Gravity application #
########################
remove recorders
eval "recorder Node -file soilDispX.out -time -node $NodeList -dof 1 disp"
eval "recorder Node -file soilDispY.out -time -node $NodeList -dof 2 disp"
eval "recorder Node -file soilDispZ.out -time -node $NodeList -dof 3 disp"
# record pore water pressure in soil elements
eval "recorder Node -file pwp.out -time -node $NodeList -dof 4 vel"
# record element stresses
eval "recorder Element -file soilstress.out -time -ele $SoilEleList material 1 stress"
# record contact element stresses
eval "recorder Element -file contactstress.out -time -ele $contactEleList material 9 stress"
# elastic behavior ---------------------------------------------- 1st run
updateMaterialStage -material 1 -stage 0
updateMaterialStage -material 9 -stage 0
set v 1.5
numberer RCM
system SparseSPD
#test NormDispIncr 1.00e-004 50 1
test NormDispIncr 1.00e-003 50 1
#test EnergyIncr 1.00e-004 50 1
algorithm ModifiedNewton
#algorithm KrylovNewton
constraints Penalty 1.e18 1.e18
#integrator Newmark $v [expr pow($v+0.5, 2)/4]
integrator Newmark 1.5 1.
analysis Transient
# preform one time step increment with a size of 5000
set ok [analyze 10 5000]
if {$ok != 0} { return $ok }
puts "First run done."
# switch material stage from elastic (gravity) to plastic ------- 2nd run
updateMaterialStage -material 1 -stage 1
updateMaterialStage -material 9 -stage 1
# perform one time step increment with a size of 1
set ok [analyze 5 1]
#set ok [analyze 1 1]
#set ok [analyze 500 0.002] #when convergence criteria =10e-004
if {$ok != 0} { return $ok }
puts "Second run done."
###################################################################################################################
# start third run with beam contact elements
###################################################################################################################
wipeAnalysis
loadConst -time 0.
##########
# Set up pile L-K damping
##########
model basic -ndm 3 -ndf 3
# set up Lysmer-Kuhlemeyer dashpots for pile
set pshearcoeff 20
set pnormcoeff 20
# add Lysmer-Kuhlemeyer damping at pile tip
uniaxialMaterial Viscous 6 $pnormcoeff 1
uniaxialMaterial Viscous 7 $pshearcoeff 1
# normal direction dashpots
node 800 0.0 -9.0 0.0
node 801 0.0 -9.0 0.0
fix 800 1 1 1
fix 801 1 0 1
element zeroLength 350 800 801 -mat 6 -dir 2
# shear direction dashpots
node 802 0.0 -9.0 0.0
node 803 0.0 -9.0 0.0
fix 802 1 1 1
fix 803 0 1 1
element zeroLength 351 802 803 -mat 7 -dir 1
wipeAnalysis
########################
# Set up pile analysis #
########################
model basic -ndm 3 -ndf 6
# pile nodes
node 1005 0.0 -9.0 0.0
node 1006 0.0 0.0 0.0
node 1007 0.0 -3.0 0.0
node 1008 0.0 -6.0 0.0
set PileNodeList {1005 1006 1007 1008}
# fix pile nodes from moving out of the symmetrical plane
#fix 1006 0 0 1 1 1 1
#fix 1007 0 0 1 1 1 1
#fix 1008 0 0 1 1 1 1
fix 1006 0 0 1 1 1 0
fix 1007 0 0 1 1 1 0
fix 1008 0 0 1 1 1 0
# fix base of pile from movement in symmetrical plane
#fix 1005 0 0 1 1 1 1
fix 1005 0 0 1 1 1 0
#fix 1005 1 0 1 1 1 1
#fix 1005 0 0 0 0 0 0
# restrain the pile from downwards vertical movement
uniaxialMaterial ENT 3 1000000
node 900 0.0 -9.0 0.0
fix 900 1 0 1 1 1 1
equalDOF 801 900 1 2 3
#equalDOF 801 900 2
element zeroLength 500 900 1005 -mat 3 -dir 2
# restrain pile in lateral direction
equalDOF 803 1005 1 2 3
#equalDOF 803 1005 1
set numSteps3rdRun 5
#mass 1006 0.0954 0.0954 0.0954 0.0954 0.0954 0.0954
mass 1007 0.1909 0.1909 0.1909 0.1909 0.1909 0.1909
mass 1008 0.1909 0.1909 0.1909 0.1909 0.1909 0.1909
#mass 1005 0.0954 0.0954 0.0954 0.0954 0.0954 0.0954
pattern Plain 101 Linear {
# load factor times reference load is the load applied to the pile
# load 1006 0. [expr -0.0954*$g/$numSteps3rdRun] 0. 0. 0. 0.
load 1007 0. [expr -0.1909*$g/$numSteps3rdRun] 0. 0. 0. 0.
load 1008 0. [expr -0.1909*$g/$numSteps3rdRun] 0. 0. 0. 0.
# load 1005 0. [expr -0.0954*$g/$numSteps3rdRun] 0. 0. 0. 0.
}
# rigidLink nodes in identical position
node 557 -0.29 0.0 0.0
node 558 0.0 0.0 0.29
node 559 0.29 0.0 0.0
node 560 -0.29 -3.0 0.0
node 561 0.0 -3.0 0.29
node 562 0.29 -3.0 0.0
node 563 -0.29 -6.0 0.0
node 564 0.0 -6.0 0.29
node 565 0.29 -6.0 0.0
node 566 -0.29 -9.0 0.0
node 567 0.0 -9.0 0.29
node 568 0.29 -9.0 0.0
node 5100 [expr 0.29*sin(45*$deg)] 0.0 [expr 0.29*sin(45*$deg)]
node 5101 [expr -0.29*sin(45*$deg)] 0.0 [expr 0.29*sin(45*$deg)]
node 5102 [expr 0.29*sin(45*$deg)] -3.0 [expr 0.29*sin(45*$deg)]
node 5103 [expr -0.29*sin(45*$deg)] -3.0 [expr 0.29*sin(45*$deg)]
node 5104 [expr 0.29*sin(45*$deg)] -6.0 [expr 0.29*sin(45*$deg)]
node 5105 [expr -0.29*sin(45*$deg)] -6.0 [expr 0.29*sin(45*$deg)]
node 5106 [expr 0.29*sin(45*$deg)] -9.0 [expr 0.29*sin(45*$deg)]
node 5107 [expr -0.29*sin(45*$deg)] -9.0 [expr 0.29*sin(45*$deg)]
# fix the link nodes in the symmetrical plane
#fix 557 0 0 1 1 1 1
#fix 559 0 0 1 1 1 1
#fix 560 0 0 1 1 1 1
#fix 562 0 0 1 1 1 1
#fix 563 0 0 1 1 1 1
#fix 565 0 0 1 1 1 1
#fix 566 0 0 1 1 1 1
#fix 568 0 0 1 1 1 1
fix 557 0 0 1 1 1 0
fix 559 0 0 1 1 1 0
fix 560 0 0 1 1 1 0
fix 562 0 0 1 1 1 0
fix 563 0 0 1 1 1 0
fix 565 0 0 1 1 1 0
fix 566 0 0 1 1 1 0
fix 568 0 0 1 1 1 0
# fix 3dof nodes to 4dof nodes at identical location
#equalDOF 557 57 1 2 3
#equalDOF 558 58 1 2 3
#equalDOF 559 59 1 2 3
#equalDOF 560 60 1 2 3
#equalDOF 561 61 1 2 3
#equalDOF 562 62 1 2 3
#equalDOF 563 63 1 2 3
#equalDOF 564 64 1 2 3
#equalDOF 565 65 1 2 3
#equalDOF 566 66 1 2 3
#equalDOF 567 67 1 2 3
#equalDOF 568 68 1 2 3
#equalDOF 5100 100 1 2 3
#equalDOF 5101 101 1 2 3
#equalDOF 5102 102 1 2 3
#equalDOF 5103 103 1 2 3
#equalDOF 5104 104 1 2 3
#equalDOF 5105 105 1 2 3
#equalDOF 5106 106 1 2 3
#equalDOF 5107 107 1 2 3
equalDOF 57 557 1 2 3
equalDOF 58 558 1 2 3
equalDOF 59 559 1 2 3
equalDOF 60 560 1 2 3
equalDOF 61 561 1 2 3
equalDOF 62 562 1 2 3
equalDOF 63 563 1 2 3
equalDOF 64 564 1 2 3
equalDOF 65 565 1 2 3
equalDOF 66 566 1 2 3
equalDOF 67 567 1 2 3
equalDOF 68 568 1 2 3
equalDOF 100 5100 1 2 3
equalDOF 101 5101 1 2 3
equalDOF 102 5102 1 2 3
equalDOF 103 5103 1 2 3
equalDOF 104 5104 1 2 3
equalDOF 105 5105 1 2 3
equalDOF 106 5106 1 2 3
equalDOF 107 5107 1 2 3
# generate pile
#geomTransf Linear 1001 0 0 -1
geomTransf Linear 1001 -1 0 0
# pile properties - half round area on line of symmetry
set I 0.0032
set A 0.1414
set E 1002320
set G 103180
set Jx 0.0064
element elasticBeamColumn 101 1007 1006 $A $E $G $Jx $I $I 1001
element elasticBeamColumn 102 1008 1007 $A $E $G $Jx $I $I 1001
element elasticBeamColumn 103 1005 1008 $A $E $G $Jx $I $I 1001
set PileElementList {101 102 103}
# set up rigid links
#geomTransf Linear 1002 -1 0 0
#geomTransf Linear 1002 0 -1 0 #also works
geomTransf Linear 1002 0 1 0
set I [expr 0.0032*1000]
element elasticBeamColumn 400 1006 557 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 401 1006 5101 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 402 1006 558 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 403 1006 5100 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 404 1006 559 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 405 1007 560 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 406 1007 5103 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 407 1007 561 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 408 1007 5102 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 409 1007 562 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 410 1008 563 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 411 1008 5105 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 412 1008 564 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 413 1008 5104 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 414 1008 565 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 415 1005 566 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 416 1005 5107 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 417 1005 567 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 418 1005 5106 $A $E $G $Jx $I $I 1002
element elasticBeamColumn 419 1005 568 $A $E $G $Jx $I $I 1002
# set up recorders
# record soil and pile displacements
eval "recorder Node -file piledispX.out -time -node $PileNodeList -dof 1 disp"
eval "recorder Node -file piledispY.out -time -node $PileNodeList -dof 2 disp"
eval "recorder Node -file piledispZ.out -time -node $PileNodeList -dof 3 disp"
# record pile reaction
eval "recorder Element -file pilereaction.out -time -ele $PileElementList globalForce"
# record accelerations
eval "recorder Node -file acceleration1.out -time -node $NodeList -dof 1 accel"
eval "recorder Node -file acceleration2.out -time -node $NodeList -dof 2 accel"
eval "recorder Node -file acceleration3.out -time -node $NodeList -dof 3 accel"
# add pile gravity (plastic) ----------- 3nd run
numberer RCM
system SparseGeneral
#test NormDispIncr 1.00e-004 50 1
test EnergyIncr 1.00e-004 50 1
algorithm KrylovNewton
constraints Penalty 1.e12 1.e12
# integrator Newmark 1.5 1.
integrator Newmark 0.5 0.25
#integrator Newmark $v [expr pow($v+0.5, 2)/4]
analysis Transient
# changed numSteps from 5 to 1
set ok [analyze $numSteps3rdRun 1]
if {$ok != 0} { return $ok }
puts "Third run done."
#############################
# Loading Earthquake
#############################
wipeAnalysis
loadConst -time 0.0
model basic -ndm 3 -ndf 6
# base input motion
pattern UniformExcitation 1 1 -accel "Series -fileTime acc_time.dat -filePath acc_value.dat -factor [expr 1.0*$g]"
# Update soil permeability properties
for {set i 10001} {[expr $i-10000]<=$NumElements} {incr i 1} {
parameter $i element [expr $i-10000] hPerm
updateParameter $i [expr 1.0000e-004/$g/$fmass]
}
for {set j 20001} {[expr $j-20000]<=$NumElements} {incr j 1} {
parameter $j element [expr $j-20000] vPerm
updateParameter $j [expr 1.0000e-004/$g/$fmass]
}
# use large damping for pushover analysis
set am 0.21542
set ak 0.00090946
set numSteps 2000
set TimeIncr 0.01
set TimeIncrMin [expr $TimeIncr/100]
set finalTime [expr $numSteps*$TimeIncr]
set gamma 0.600000
numberer RCM
system SparseGeneral
#test EnergyIncr 1.00e-004 10 1
test NormDispIncr 1.00e-004 50 1
#algorithm KrylovNewton
algorithm ModifiedNewton
constraints Penalty 1.e12 1.e12
rayleigh $am 0.0 $ak 0.0
# can use either integrator
integrator Newmark $gamma [expr pow($gamma+0.5, 2)/4]
#integrator Newmark 0.5 0.25
analysis VariableTransient
#analysis Transient
set ok 0
set time [getTime]
while {$time<$finalTime && $ok==0} {
set ok [analyze 1 $TimeIncr $TimeIncrMin $TimeIncr 15]
set time [getTime]
puts "$time"
}
#analyze $numSteps $TimeIncr [expr $TimeIncr/100] $TimeIncr 15
#analyze $numSteps $TimeIncr
puts "Analysis Finished!!"
wipe
Convergence followed by divergence in the same analysis
Moderator: Moderators