hello every one , I am modeling a buried pipeline in opensees and to perform dynamic analysis on it.I have checked my source code by partially running it and it works fine but when i do analysis, it does not run and start giving following errors, Does anybody have an idea why this is happening.Please Help
WARNING SuperLU::solve(void)-Error 5 returned in factorization dgstrf()
WARNING NewtonRaphson::solveCurrentStep()-the LinearSysOfEqn failed in solve()
DirectIntegrationAnalysis::analyze()-the algorithm failed at time 0.00976
OpemSees > analyze failed, returned: -3 error flag
Here is my source Code
# Dynamic analysis of a single pipeline
# As Beam on nonlinear Wnkler Foundation
# Lateral pipe response is described by p-y springs.
# Basic units: KN, m
#
########################
wipe
set pi 3.14159265358979
set gravity 9.81
########################
#pipe geometry and mesh
########################
# Length of embedded pipe (below ground surface) (m)
set L1 25.0
#pipe diameter(m)
set D 1.0
#pipe thickness(m)
set t 0.04
#pipe area(m^2)
set A [expr $pi / 4 * (pow($D, 2)- pow($D-2*$t, 2))]
puts "Area is $A"
#pipe stiffness(kN/m^2)
#set E 2.1e8
set E 6.92e7
set Emass [expr $E * 1000]
#Moment of inertia(m^4)
set I [expr $pi / 64 * (pow($D, 4) - pow($D-2*$t, 4))]
puts "Inertia moment is $I"
puts "[expr $E*$I]"
#number of pipe elements
set nElePipe 46
#pipe element length
set eleSize [expr ($L1) / $nElePipe]
#number of total pipe nodes
set nNodePipe [expr 1 + $nElePipe]
#---RAYLEIGH DAMPING PARAMETERS
# damping ratio
set damp 0.02
# lower frequency
set omega1 [expr 2*$pi*1.5]
# upper frequency
set omega2 [expr 2*$pi*20]
puts "done"
# damping coefficients
set alpha1 [expr 2*$damp*$omega1*$omega2/($omega1 + $omega2)]
set beta1 [expr 2*$damp/($omega1 + $omega2)]
puts "damping coefficients: alpha = $alpha1; beta = $beta1"
########################
#Soil properties
########################
#soil unit weigth (kN/m^3)
set gamma 15.8
#soil friction angle (degrees)
set phi 40.0
# select pult definition method for p-y curves
# API (default) --> 1
# Brinch Hansen --> 2
set puSwitch 1
# variation in coefficent of subgrade reaction with depth for p-y curves
# API linear variation (default) --> 1
# modified API parabolic variation --> 2
set kSwitch 1
# effect of ground water on subgrade reaction modulus for p-y curves
# above gwt --> 1
# below gwt --> 2
set gwtSwitch 1
#shear wave velocity (m/s) : avg. value at mid. depth(11m)
set Vs 290
########################
#Node
########################
#######################
#create pile nodes
#######################
model BasicBuilder -ndm 2 -ndf 3
for {set i 1} {$i <= $nNodePipe} {incr i} {
set yCoord [expr $eleSize * ($i - 1)]
node $i 0.00 $yCoord
}
puts "Finished creating all pipe nodes..."
#create fixity at pipe base
fix 1 0 1 0
#fix 1 1 1 0
# nodal masses
#mass 21 300 1.e-9 0.; # node#, Mx My Mz,
mass 20 300 1.e-9 0.; # node#, Mx My Mz,
mass 21 300 1.e-9 0.; # node#, Mx My Mz,
mass 22 300 1.e-9 0.; # node#, Mx My Mz,
puts "Finished creating all pipe node fixities and mass..."
##########################
#create spring nodes
##########################
model BasicBuilder -ndm 2 -ndf 3
set count 0
for {set i 1} {$i<= $nNodePipe} {incr i} {
set yCoord [expr $eleSize * ($i - 1)]
if {$yCoord <= $L1-$eleSize} {
node [expr $i + 100] 0.00 $yCoord # slave node
node [expr $i + 200] 0.00 $yCoord # fixed node
set count [expr $count + 1]
}
}
puts "Finished creating all spring nodes..."
#number of Embedded nodes
set nNodeEmbed $count
puts "embeded number: $nNodeEmbed"
#spring nodes fixities
for {set i 1} {$i <= $nNodeEmbed} {incr i} {
fix [expr $i + 100] 1 0 1 #slave node
fix [expr $i + 200] 1 1 1 #fixed node
}
#create coordinate-transformation object
geomTransf Linear 1
puts "Finished creating all spring node fixities..."
########################
#define equal dof between pile and spring nodes
########################
for {set i 1} {$i <= $nNodeEmbed} {incr i} {
equalDOF $i [expr $i + 100] 1
}
puts "Finished creating all equal degrees of freedom..."
########################
#Create spring material
########################
#specify necessary procedures
#source get_pyParamji_final.tcl
#p-y spring matrial
#set new py material parameters
### set parameters for MT equation (unit is N, cm)
# atmosphere pressure
set Pa 10.13
# soil unit weight
set gamma2 [expr $gamma/1000]
# passive pressure coeff.
set Kp [expr pow(tan([expr (45+$phi/2)*$pi/180 ]),2)]
# static pressure coeff.
set K0 [expr 1 - sin([expr $phi*$pi/180])]
for {set i 1} {$i <= $nNodeEmbed} {incr i} {
#depth of current py node
set pyDepth [expr ($L1 - $eleSize * ($i - 1))]
set pult [expr 13.3*$D*100*$Kp*$gamma2*pow($pyDepth*100,1.02)]
set pult [expr $pult*0.1]
set pult_U [expr $pult*$eleSize]
puts "at node $i pult = [expr $pult]"
# Determine maximum subgrade reaction modulus (kmax) based on Mylonakis study
set Ep [expr 6.92E7*0.284]
set stressratio [expr $pyDepth * $gamma / 101.3]
set Vs [expr 350 * pow($stressratio, 0.282)]
set G [expr $gamma / 9.81 * pow($Vs, 2)]
set Es [expr 2 * $G * (1 + 0.3)]
set stiffratio [expr $Ep / $Es]
set kmax [expr 1.67 * pow($stiffratio, -0.053) * $Es]
# set kmax [expr 6.918 * pow($stiffratio, -0.137) * $Es]
puts "at node $i kmax = $kmax"
set kmax_U [expr $kmax * 0.5]
# Determine yielding strain based on Darendeli thesis
set mean [expr $pyDepth * $gamma * (1+2*$K0)/3]
set gamma_r [expr 0.0475 * pow($mean, 0.239)]
set gamma_y [expr $gamma_r * pow([expr 1/0.99 -1], 1/0.895)]
set gamma_y [expr $gamma_y / 100]
puts "at node $i gamma_y = $gamma_y"
# Determine yielding disp. based on elastic FEA results for round shape pipe
set y_yield [expr $gamma_y * $D / 0.4]
# Determine yielding resistance
set p_yield_U [expr $kmax * $y_yield * $eleSize]
puts "at node $i p_yield = $p_yield_U"
# Determine plastic stiffness term, H
set H [expr 0.015 * pow([expr $mean / 101.3], -0.564) * $kmax_U]
# set H [expr 0.1 / 3.5 * pow([expr $mean / 101.3], -0.564) * $kmax_U]
puts "at node $i H = [expr $H]"
# tag pult pyield ke H
uniaxialMaterial PySimple1 [expr $i + 100 ] 1 [expr $pult_U] $p_yield_U $kmax_U $H
# uniaxialMaterial PySimple3 [expr $i + 100] [expr $pult*0.5] [expr $p_yield*0.5] [expr $kmax] $H
}
########################
#create zero-length elements for spring
########################
for {set i 1} {$i <= $nNodeEmbed} {incr i} {
element zeroLength [expr $i + 1000] [expr $i + 100] [expr $i + 200] -mat [expr $i +100] -dir 1
}
puts "Finished creating all zerolength elements..."
#set f [eigen 2]
#puts "f is $f"
set f [expr [eigen fullGenLapack 1]**0.5]
set f [expr $f / 2 / $pi]
puts "$f"
puts "Finished creating all p-y springs..."
########################
#create pile elements
########################
for {set i 1} {$i <= 46} {incr i} {
#eletag #iNode #jNode #Area #E #Iz #Transtag
element elasticBeamColumn $i $i [expr $i + 1] $A $E $I 1
}
puts "Finished creating all pipe elements..."
########################
#create recorders
########################
#set timeStep 0.1
recorder Node -file pipeDisp_Large_0.1g_gMotion_d.out -time -nodeRange 1 $nNodePipe -dof 1 disp
recorder Node -file pipe_acc_Large_0.1g_gMotion_d.out -time -node 1 $nNodePipe -dof 1 accel
################################
#Dynamic Ground-Motion Analysis
################################
#timeSeries Path 2 -dt 0.00976 -filePath sine_cen_0.3g.tcl -factor $gravity; # define acceleration vector from file (dt=0.00976 is associated with the input file gm)
for {set i 1} {$i<= $nNodeEmbed} {incr i} {
timeSeries Path [expr $i+100] -dt 0.00976 -filePath accel/acc_node$i.txt -factor 1
timeSeries Path [expr $i+200] -dt 0.00976 -filePath velocity/vel_node$i.txt -factor 1
timeSeries Path [expr $i+300] -dt 0.00976 -filePath disp/disp_node$i.txt -factor 1
}
pattern MultipleSupport 101 {
for {set i 1} {$i<= $nNodeEmbed} {incr i} {
groundMotion $i Plain -accel [expr $i+100] -vel [expr $i+200] -disp [expr $i+300]
imposedMotion [expr $i+200] 1 $i #node, dof, GM tag
}
}
# set damping based on first eigen mode
set freq 1.5
#[expr pow($det, 0.5)]
set dampRatio 0.02
set beta [expr 2*$dampRatio/$freq]
#rayleigh 0. 0. 0. [expr 2*$dampRatio/$freq]
#rayleigh 0. 0. 0. $beta1
#rayleigh 0.05 0.0003 0.0 0.0
rayleigh $alpha1 $beta1 0.0 0.0
# create the analysis
#wipeAnalysis; # clear previously-define analysis parameters
#constraints Plain;
#constraints Transformation; # how it handles boundary conditions
constraints Penalty 1.e18 1.e18
#test NormDispIncr 1e-12 20
test NormDispIncr 1e-9 20
numberer RCM; # renumber dof's to minimize band-width (optimization), if you want to
#system UmfPack
system BandGeneral; # how to store and solve the system of equations in the analysis
algorithm ModifiedNewton # use Linear algorithm for linear analysis
integrator Newmark 0.5 0.25 ; # determine the next time step for an analysis
#analysis Transient; # define type of analysis: time-dependent
analysis VariableTransient; # define type of analysis: time-dependent
set startT [clock seconds]
puts "Starting Load Application..."
#analyze 2049 0.01952 # apply 2049 0.01952-sec time steps in analysis
#analyze 2049 0.01952 [expr 0.01952/64.0] 0.01952 10; # no. of time step, time step incr., dtmin, dtmax, no. of interation
analyze 1560 0.00976;
set endT [clock seconds]
puts "Load Application finished..."
puts "Loading Analysis execution time: [expr $endT-$startT] seconds."
Pipeline Modeling
Moderators: silvia, selimgunay, Moderators
Re: Pipeline Modeling
1. I think you can change nodes fixity like that
#spring nodes fixities
for {set i 1} {$i <= $nNodeEmbed} {incr i} {
fix [expr $i + 100] 0 1 1
fix [expr $i + 200] 1 1 1
}
2. You can increase tolerance on the norm of the displacement residual was 1e3
test NormDispIncr 1e3 20
I've changed your code, and it can be run. But, I don't know why do you put y50=$p_yield_U, cd=$kmax_U, and c=$H in your model.
#spring nodes fixities
for {set i 1} {$i <= $nNodeEmbed} {incr i} {
fix [expr $i + 100] 0 1 1
fix [expr $i + 200] 1 1 1
}
2. You can increase tolerance on the norm of the displacement residual was 1e3
test NormDispIncr 1e3 20
I've changed your code, and it can be run. But, I don't know why do you put y50=$p_yield_U, cd=$kmax_U, and c=$H in your model.