Pipeline Modeling

[billu]

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."

[Baodkt]

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.