Hi Dear all,
I am building a 2D soil model with 9_4_quadup element and PressureDependYield02 material. I have a problem. During the analysis, these errors come up
WARNING: AcceleratedNewton::solveCurrentStep()the LinearSysofEqn failed in solve ()
NewtonRaphson:: solveCurrentStep() -the Convergence Test object failed in test()
DirectIntegration Analysis : : analyze() - the Algorithm failed at time 50000
OpenSees > analyze failed, returned: -3 error flag
It seems the analysis failed at the first step of elastic gravity analysis stage. I am wondering if it is caused by the inadequate constraints in the model. But I have checked my tcl file many times and cannot find the bug. I have posted my model as below, could you please give me a hand to fix it? I appreciate in advance.
wipe
#-----------------------------------------------------------------------------------------
# 1. DEFINE SOIL AND MESH GEOMETRY
#-----------------------------------------------------------------------------------------
#---SOIL GEOMETRY
# thicknesses of soil profile (m)
set soilThick 5.0
# number of soil layers
set numLayers 2
# layer thicknesses
set layerThick(2) 3.0
set layerThick(1) 2.0
# depth of water table
set waterTable 0.0
# define layer boundaries
set layerBound(1) $layerThick(1)
for {set i 2} {$i <= $numLayers} {incr i 1} {
set layerBound($i) [expr $layerBound([expr $i-1])+$layerThick($i)]
}
#---MESH GEOMETRY
# number of elements in horizontal direction
set nElemX(1) 5
set nElemX(2) 6
set nElemX(3) 5
# total number of elements in vertical direction
set numElemX 16
# number of nodes in horizontal direction
set nNodeX(1) 11
set nNodeX(2) 23
set nNodeX(3) 33
set numNodeX [expr 2*$numElemX+1]
# horizontal element size (m)
set sElemX(1) 0.7
set sElemX(2) 0.5
set sElemX(3) 0.7
#number of elements in vertical direction for each layer
set nElemY(2) 15
set nElemY(1) 5
# total number of elements in vertical direction
set nElemT 20
# vertical element size in each layer
for {set i 1} {$i <=$numLayers} {incr i 1} {
set sElemY($i) [expr $layerThick($i)/$nElemY($i)]
puts "size: $sElemY($i)"
}
# number of nodes in vertical direction
set nNodeY [expr 2*$nElemT+1]
# total number of nodes
set nNodeT [expr $numNodeX*$nNodeY]
#-----------------------------------------------------------------------------------------
# 2. CREATE PORE PRESSURE NODES AND FIXITIES
#-----------------------------------------------------------------------------------------
model BasicBuilder -ndm 2 -ndf 3
set ppNodesInfo [open ppNodesInfo.dat w]
set count 1
set layerNodeCount 0
# loop over soil layers
for {set k 1} {$k <= $numLayers} {incr k 1} {
# loop in horizontal direction
for {set i 1} {$i <= $numNodeX} {incr i 2} {
# loop in vertical direction
if {$k == 1} {
set bump 1
} else {
set bump 0
}
for {set j 1} {$j <= [expr 2*$nElemY($k)+$bump]} {incr j 2} {
if {$i <= $nNodeX(1)} {
set xCoord [expr ($i-1)*$sElemX(1)/2]
} elseif {$i > $nNodeX(1) && $i <= $nNodeX(2)} {
set xCoord [expr 3.5+($i-$nNodeX(1))*$sElemX(2)/2]
} else {
set xCoord [expr 6.5+($i-$nNodeX(2))*$sElemX(3)/2]
}
set yctr [expr $j + $layerNodeCount]
if {$yctr <= [expr 2*$nElemY(1)+1]} {
set yCoord [expr ($yctr-1)*$sElemY($k)/2]
} else {
set yCoord [expr $layerThick(1)+($yctr-11)*$sElemY($k)/2]
}
set nodeNum [expr $i + ($yctr-1)*$numNodeX]
node $nodeNum $xCoord $yCoord
# output nodal information to data file
puts $ppNodesInfo "$nodeNum $xCoord $yCoord"
# designate nodes above water table
set waterHeight [expr $soilThick-$waterTable]
if {$yCoord>=$waterHeight} {
set dryNode($count) $nodeNum
set count [expr $count+1]
}
}
}
set layerNodeCount [expr $yctr + 1]
}
close $ppNodesInfo
puts "Finished creating all -ndf 3 nodes..."
# define fixities for pore pressure nodes above water table
for {set i 1} {$i < $count} {incr i 1} {
fix $dryNode($i) 0 0 1
}
# define fixities for pore pressure nodes at base of soil column
fix 1 0 1 0
fix 3 0 1 0
fix 5 0 1 0
fix 7 0 1 0
fix 9 0 1 0
fix 11 0 1 0
fix 13 0 1 0
fix 15 0 1 0
fix 17 0 1 0
fix 19 0 1 0
fix 21 0 1 0
fix 23 0 1 0
fix 25 0 1 0
fix 27 0 1 0
fix 29 0 1 0
fix 31 0 1 0
fix 33 0 1 0
puts "Finished creating all -ndf 3 boundary conditions..."
# define equal degrees of freedom for pore pressure nodes
for {set i 1} {$i <= [expr $numNodeX-2]} {incr i 2} {
for {set j 1} {$j <= $nNodeY} {incr j 2} {
equalDOF [expr $i+($j-1)*33] [expr $i+2+($j-1)*33] 1 2
}
}
puts "Finished creating equalDOF for pore pressure nodes..."
#-----------------------------------------------------------------------------------------
# 3. CREATE INTERIOR NODES AND FIXITIES
#-----------------------------------------------------------------------------------------
model BasicBuilder -ndm 2 -ndf 2
# central column of nodes
set layerNodeCount 0
# loop over soil layers
for {set k 1} {$k <= $numLayers} {incr k 1} {
# loop in horizontal direction
for {set i 2} {$i <= $numNodeX} {incr i 2} {
# loop in vertical direction
if {$k == 1} {
set bump 1
} else {
set bump 0
}
for {set j 1} {$j <= [expr 2*$nElemY($k)+$bump]} {incr j 1} {
if {$i <= $nNodeX(1)} {
set xCoord [expr ($i-1)*$sElemX(1)/2]
} elseif {$i > $nNodeX(1) && $i <= $nNodeX(2)} {
set xCoord [expr 3.5+($i-$nNodeX(1))*$sElemX(2)/2]
} else {
set xCoord [expr 6.5+($i-$nNodeX(2))*$sElemX(3)/2]
}
set yctr [expr $j + $layerNodeCount]
if {$yctr <= [expr 2*$nElemY(1)+1]} {
set yCoord [expr ($yctr-1)*$sElemY($k)/2]
} else {
set yCoord [expr $layerThick(1)+($yctr-11)*$sElemY($k)/2]
}
set nodeNum [expr $i + ($yctr-1)*$numNodeX]
node $nodeNum $xCoord $yCoord
}
}
set layerNodeCount $yctr
}
puts "Finished creating central column nodes..."
# interior nodes on the element edges
set layerNodeCount 0
# loop over layers
for {set k 1} {$k <= $numLayers} {incr k 1} {
# loop in horizontal direction
for {set i 1} {$i <= $numNodeX} {incr i 2} {
# loop in vertical direction
for {set j 1} {$j <= $nElemY($k)} {incr j 1} {
if {$i <= $nNodeX(1)} {
set xCoord [expr ($i-1)*$sElemX(1)/2]
} elseif {$i > $nNodeX(1) && $i <= $nNodeX(2)} {
set xCoord [expr 3.5+($i-$nNodeX(1))*$sElemX(2)/2]
} else {
set xCoord [expr 6.5+($i-$nNodeX(2))*$sElemX(3)/2]
}
set yctr [expr $j + $layerNodeCount]
if {$yctr <= $nElemY(1)} {
set yCoord [expr $sElemY($k)*($yctr-0.5)]
} else {
set yCoord [expr $layerThick(1)+$sElemY($k)*($yctr-5.5)]
}
set nodeNum [expr $i + (2*$yctr-1)*$numNodeX]
node $nodeNum $xCoord $yCoord
}
}
set layerNodeCount $yctr
}
puts "Finished creating all -ndf 2 nodes..."
# define fixities for interior nodes at base of soil column
fix 2 0 1
fix 4 0 1
fix 6 0 1
fix 8 0 1
fix 10 0 1
fix 12 0 1
fix 14 0 1
fix 16 0 1
fix 18 0 1
fix 20 0 1
fix 22 0 1
fix 24 0 1
fix 26 0 1
fix 28 0 1
fix 30 0 1
fix 32 0 1
puts "Finished creating all -ndf 2 boundary conditions..."
# define equal degrees of freedom which have not yet been defined
for {set i 1} {$i <= [expr $numNodeX-2]} {incr i 2} {
for {set j 1} {$j <= $nNodeY} {incr j 1} {
equalDOF [expr $i+($j-1)*33] [expr $i+1+($j-1)*33] 1 2
}
}
puts "Finished creating equalDOF constraints..."
#-----------------------------------------------------------------------------------------
## 4. CREATE SOIL MATERIALS
#-----------------------------------------------------------------------------------------
# define grade of slope (%)
set grade 0.0
set slope [expr atan($grade/100.0)]
set g -9.81
nDMaterial PressureDependMultiYield02 2 2 1.9 1.0e5 2.33e5 33.5 0.1 \
101.0 0.5 25.5 0.045 0.15 0.06 \
0.15 20 5.0 3.0 1.0 \
0.0 0.7 0.9 0.02 0.7 101.0
set thick(2) 0.5
set xWgt(2) [expr $g*sin($slope)]
set yWgt(2) [expr $g*cos($slope)]
set uBulk(2) 2.2e6
set hPerm(2) [expr 6.6000e-005/9.81]
set vPerm(2) [expr 6.6000e-005/9.81]
nDMaterial PressureDependMultiYield02 1 2 2.1 1.3e5 2.6e5 36.5 0.1 \
101.0 0.5 26 0.013 0.0 0.3 \
0.0 20 5.0 3.0 1.0 \
0.0 0.55 0.9 0.02 0.7 101.0
set thick(1) 0.5
set xWgt(1) [expr $g*sin($slope)]
set yWgt(1) [expr $g*cos($slope)]
set uBulk(1) 2.2e6
set hPerm(1) [expr 6.6000e-005/9.81]
set vPerm(1) [expr 6.6000e-005/9.81]
puts "Finished creating all soil materials..."
#-----------------------------------------------------------------------------------------
# 5. CREATE SOIL ELEMENTS
#-----------------------------------------------------------------------------------------
set ElemInfo [open ElemInfo.dat w]
for {set j 1} {$j <= $nElemT} {incr j 1} {
for {set i 1} {$i <= $numElemX} {incr i 1} {
set nI [expr (2*$i-1)+($j-1)*66]
set nJ [expr $nI + 2]
set nK [expr $nI + 68]
set nL [expr $nI + 66]
set nM [expr $nI + 1]
set nN [expr $nI + 35]
set nP [expr $nI + 67]
set nQ [expr $nI + 33]
set nR [expr $nI + 34]
set lowerBound 0.0
for {set k 1} {$k <= $numLayers} {incr k 1} {
if {[expr $j*$sElemY($k)] <= $layerBound($k) && [expr $j*$sElemY($k)] > $lowerBound} {
# permeabilities are initially set at 100/9.81 m/s for gravity analysis, values are updated post-gravity
element 9_4_QuadUP [expr $i+($j-1)*$numElemX] $nI $nJ $nK $nL $nM $nN $nP $nQ $nR \
$thick($k) $k $uBulk($k) 1.0 [expr 1.0000e+002/9.81] [expr 1.0000e+002/9.81] $xWgt($k) $yWgt($k)
puts $ElemInfo "[expr $i+($j-1)*$numElemX] $nI $nJ $nK $nL $nM $nN $nP $nQ $nR"
}
set lowerBound $layerBound($k)
}
}
}
close $ElemInfo
puts "Finished creating all soil elements..."
#----------------------------------------------------------------------------------------
# 7. CREATE GRAVITY RECORDERS
#-----------------------------------------------------------------------------------------
# create list for pore pressure nodes
set nodeList3 {}
set channel [open "ppNodesInfo.dat" r]
set count 0;
foreach line [split [read -nonewline $channel] \n] {
set count [expr $count+1];
set lineData($count) $line
set nodeNumber [lindex $lineData($count) 0]
lappend nodeList3 $nodeNumber
}
close $channel
# record nodal displacment, acceleration, and porepressure
eval "recorder Node -file Gdisplacement.out -time -node $nodeList3 -dof 1 2 disp"
eval "recorder Node -file Gacceleration.out -time -node $nodeList3 -dof 1 2 accel"
eval "recorder Node -file GporePressure.out -time -node $nodeList3 -dof 3 vel"
# record elemental stress and strain (files are names to reflect GiD gp numbering)
recorder Element -file Gstress1.out -time -eleRange 1 $nElemT material 1 stress
recorder Element -file Gstress2.out -time -eleRange 1 $nElemT material 2 stress
recorder Element -file Gstress3.out -time -eleRange 1 $nElemT material 3 stress
recorder Element -file Gstress4.out -time -eleRange 1 $nElemT material 4 stress
recorder Element -file Gstress9.out -time -eleRange 1 $nElemT material 9 stress
recorder Element -file Gstrain1.out -time -eleRange 1 $nElemT material 1 strain
recorder Element -file Gstrain2.out -time -eleRange 1 $nElemT material 2 strain
recorder Element -file Gstrain3.out -time -eleRange 1 $nElemT material 3 strain
recorder Element -file Gstrain4.out -time -eleRange 1 $nElemT material 4 strain
recorder Element -file Gstrain9.out -time -eleRange 1 $nElemT material 9 strain
puts "Finished creating gravity recorders..."
#-----------------------------------------------------------------------------------------
# 9. DEFINE ANALYSIS PARAMETERS
#-----------------------------------------------------------------------------------------
#---GROUND MOTION PARAMETERS
# time step in ground motion record
set motionDT 0.005
# number of steps in ground motion record
set motionSteps 5200
#---RAYLEIGH DAMPING PARAMETERS
set pi 3.141592654
# damping ratio
set damp 0.02
# lower frequency
set omega1 [expr 2*$pi*0.2]
# upper frequency
set omega2 [expr 2*$pi*20]
# damping coefficients
set a0 [expr 2*$damp*$omega1*$omega2/($omega1 + $omega2)]
set a1 [expr 2*$damp/($omega1 + $omega2)]
puts "damping coefficients: a_0 = $a0; a_1 = $a1"
#---DETERMINE STABLE ANALYSIS TIME STEP USING CFL CONDITION
# maximum shear wave velocity (m/s)
set vsMax 250.0
# duration of ground motion (s)
set duration [expr $motionDT*$motionSteps]
# minimum element size
set minSize $sElemY(1)
for {set i 2} {$i <= $numLayers} {incr i 1} {
if {$sElemY($i) < $minSize} {
set minSize $sElemY($i)
}
}
# trial analysis time step
set kTrial [expr $minSize/(pow($vsMax,0.5))]
# define time step and number of steps for analysis
if { $motionDT <= $kTrial } {
set nSteps $motionSteps
set dT $motionDT
} else {
set nSteps [expr int(floor($duration/$kTrial)+1)]
set dT [expr $duration/$nSteps]
}
puts "number of steps in analysis: $nSteps"
puts "analysis time step: $dT"
#---ANALYSIS PARAMETERS
# Newmark parameters
set gamma 0.5
set beta 0.25
#-----------------------------------------------------------------------------------------
# 10. GRAVITY ANALYSIS
#-----------------------------------------------------------------------------------------
# update materials to ensure elastic behavior
updateMaterialStage -material 1 -stage 0
updateMaterialStage -material 2 -stage 0
constraints Penalty 1.e14 1.e14
test NormDispIncr 1e-4 100 1
algorithm KrylovNewton
numberer RCM
system ProfileSPD
integrator Newmark $gamma $beta
analysis Transient
set startT [clock seconds]
set ok [analyze 5 5.00e+004]
if {$ok != 0} { return $ok }
puts "Finished with elastic gravity analysis..."
# fix base before dynamic run
model basic -ndm 2 -ndf 3
fixY [expr 0.000000e+000] 1 1 0 -tol 1e-4
# update materials to consider plastic behavior
updateMaterialStage -material 1 -stage 1
updateMaterialStage -material 2 -stage 1
# plastic gravity loading
set ok [analyze 5 5.00e+004]
if {$ok != 0} { return $ok }
puts "Finished with plastic gravity analysis..."
#-----------------------------------------------------------------------------------------
# 11. UPDATE ELEMENT PERMEABILITY VALUES FOR POST-GRAVITY ANALYSIS
#-----------------------------------------------------------------------------------------
# choose base number for parameter IDs which is higer than other tags used in analysis
set ctr 10000.0
# loop over elements to define parameter IDs
for {set i 1} {$i<=$nElemT} {incr i 1} {
parameter [expr int($ctr+1.0)] element $i vPerm
parameter [expr int($ctr+2.0)] element $i hPerm
set ctr [expr $ctr+2.0]
}
# update permeability parameters for each element using parameter IDs
set ctr 10000.0
for {set j 1} {$j <= $nElemT} {incr j 1} {
set lowerBound 0.0
for {set i 1} {$i <= $numLayers} {incr i 1} {
if {[expr $j*$sElemY($i)] <= $layerBound($i) && [expr $j*$sElemY($i)] > $lowerBound} {
updateParameter [expr int($ctr+1.0)] $vPerm($i)
updateParameter [expr int($ctr+2.0)] $hPerm($i)
}
set lowerBound $layerBound($i)
}
set ctr [expr $ctr+2.0]
}
puts "Finished updating permeabilities for dynamic analysis..."
#-----------------------------------------------------------------------------------------
# 12. CREATE POST-GRAVITY RECORDERS
#-----------------------------------------------------------------------------------------
# reset time and analysis
setTime 0.0
wipeAnalysis
remove recorders
# recorder time step
set recDT [expr 2*$motionDT]
# record nodal displacment, acceleration, and porepressure
eval "recorder Node -file displacement.out -time -dT $recDT -node $nodeList3 -dof 1 2 disp"
eval "recorder Node -file acceleration.out -time -dT $recDT -node $nodeList3 -dof 1 2 accel"
eval "recorder Node -file porePressure.out -time -dT $recDT -node $nodeList3 -dof 3 vel"
# record elemental stress and strain (files are names to reflect GiD gp numbering)
recorder Element -file stress1.out -time -dT $recDT -eleRange 1 $nElemT material 1 stress
recorder Element -file stress2.out -time -dT $recDT -eleRange 1 $nElemT material 2 stress
recorder Element -file stress3.out -time -dT $recDT -eleRange 1 $nElemT material 3 stress
recorder Element -file stress4.out -time -dT $recDT -eleRange 1 $nElemT material 4 stress
recorder Element -file stress9.out -time -dT $recDT -eleRange 1 $nElemT material 9 stress
recorder Element -file strain1.out -time -dT $recDT -eleRange 1 $nElemT material 1 strain
recorder Element -file strain2.out -time -dT $recDT -eleRange 1 $nElemT material 2 strain
recorder Element -file strain3.out -time -dT $recDT -eleRange 1 $nElemT material 3 strain
recorder Element -file strain4.out -time -dT $recDT -eleRange 1 $nElemT material 4 strain
recorder Element -file strain9.out -time -dT $recDT -eleRange 1 $nElemT material 9 strain
puts "Finished creating all recorders..."
#-----------------------------------------------------------------------------------------
# 13. DYNAMIC ANALYSIS
#-----------------------------------------------------------------------------------------
model BasicBuilder -ndm 2 -ndf 3
loadConst -time 0.0
# base input motion
pattern UniformExcitation 1 1 -accel "Series -fileTime acc_time.dat -filePath acc_value.dat -factor [expr 1.00000e+000*9.81]"
puts "Dynamic loading created..."
constraints Penalty 1.e16 1.e16
test NormDispIncr 1.0e-3 35 1
algorithm KrylovNewton
numberer RCM
system ProfileSPD
integrator Newmark $gamma $beta
rayleigh $a0 $a1 0.0 0.0
analysis Transient
# perform analysis with timestep reduction loop
set ok [analyze $nSteps $dT]
# if analysis fails, reduce timestep and continue with analysis
if {$ok != 0} {
puts "did not converge, reducing time step"
set curTime [getTime]
set mTime $curTime
puts "curTime: $curTime"
set curStep [expr $curTime/$dT]
puts "curStep: $curStep"
set rStep [expr ($nSteps-$curStep)*2.0]
set remStep [expr int(($nSteps-$curStep)*2.0)]
puts "remStep: $remStep"
set dT [expr $dT/2.0]
puts "dT: $dT"
set ok [analyze $remStep $dT]
# if analysis fails again, reduce timestep and continue with analysis
if {$ok != 0} {
puts "did not converge, reducing time step"
set curTime [getTime]
puts "curTime: $curTime"
set curStep [expr ($curTime-$mTime)/$dT]
puts "curStep: $curStep"
set remStep [expr int(($rStep-$curStep)*2.0)]
puts "remStep: $remStep"
set dT [expr $dT/2.0]
puts "dT: $dT"
analyze $remStep $dT
}
}
set endT [clock seconds]
puts "Finished with dynamic analysis..."
puts "Analysis execution time: [expr $endT-$startT] seconds"
wipe
-3 error flag
Moderators: silvia, selimgunay, Moderators