modelling interaction soil and pipeline with Beam contact 2D

mahsashamsaei · Post by **mahsashamsaei** » Wed Sep 13, 2017 10:07 pm

hi everyone
I modeled the soil and PipeLine as below, But the model is not analyzed and no output is recorded,I appreciate if anyone will analyze the model and understand the cause.
soil is a laminar box.please chek that boundry condition is tru. i used nonlinear beam column as pipeline and beam contact2D for interface element. The axial strain is the desired output.
My model is divided into two sections in the vertical direction. As sheetwall example in opensees site. and EXTRad is out diameter and width is 2*outdiameter.
When the modal analysis is carried out to determine the soil period, dont identify the contact elements.That's why this part of the modeler. for this reason remove this part from modell.
As the last question: Is the stress strain output at the section , the same strain longitudinal? or perpendicular to the strain? What is the record for getting such an output?

my Gmail address:mahsa.shamsaee@gmail.com

#############################
# #
# Modelling soil and pipe #
# #
# By #
# #
# Mahsa Shamsaei #
# #
# 1396/05/13 #
# #
#####################################

# site reponse analysis of a soil
# basic units are SI units
# thicknesses of soil profile (m)

wipe all

set parallel 0; # 0 -- sequential run (e.g., on a PC)
# 1 -- parallel run
switch $parallel {
0 {
set solver "SparseSYM"
# set solver "ProfileSPD"
# set solver "SparseGEN"
}
1 {
set solver "Mumps"
# set solver "Mumps -ICNTL14 50"
}
}

#set numberer "Plain"
set numberer "RCM"
file mkdir static
# Unit system in this file: SI

#---define UNITS ---------------------------------------------

set g 9.81; # gravitational acceleration (also good for US/EN units since SI is actually used)
set nonlinear 1

#-------------------------------------------------------------------------------------------
# 1. DEFINE ANALYSIS PARAMETERS
#-------------------------------------------------------------------------------------------
#---SOIL GEOMETRY

set extRad 0.10955
set Width [expr $extRad*2]
set V1 99.0
set V2 2.0
set V [expr $V1+$V2+$Width]
set H 200.0
#---MESH GEOMETRY
set DeltaX 2.0
set DeltaY 2.0

set NeleY1 [expr int($V1/$DeltaY)] ;# Number of element in Y dir.
set NeleY2 [expr int($V2/$DeltaY)] ;# Number of element in Y dir.
set NeleX [expr int($H/$DeltaX)] ;# Number of element in X dir.

set NNodeX [expr $NeleX+1] ;# Number of nodes in X dir.
set NNodeY1 [expr $NeleY1+1] ;# Number of nodes in Y dir.
set NNodeY2 [expr $NeleY2+1] ;# Number of nodes in Y dir.
set NNodeY [expr $NNodeY1+$NNodeY2] ;# Total Number of nodes in Y dir.

#---MATERIAL PROPERTIES
# soil mass density (Kg/m^3)
set rho 1.9
# poisson's ratio of soil
set nu 0.31
# soil shear modulus (kPa)
set refShearModul 7.5e4
# soil elastic modulus (kPa)
set Es [expr 2*$refShearModul*(1+$nu)]
# soil shear wave velocity (m/s)
set ss [expr $Es/(2*(1+$nu)*$rho)]
set Vs [expr pow($ss,0.5)]
set pp [expr ((1-$nu)*$Es)/((1-2*$nu)*(1+$nu)*$rho)]
set Vp [expr pow($pp,0.5)]
puts "Vs=$Vs Vp=$Vp"
# soil bulk modulus (kPa)
set refBulkModul 2.0e5
# soil friction angle
set frictionAng 33
set pi 3.141592654
set phi [expr ($frictionAng*$pi)/180] ;#Rad
# peak shear strain
set peakShearStra 0.1
# reference pressure
set refPress 80
# pressure dependency coefficient
set pressDependCoe 0.5
# Phase transformationangle
set PTAng 27
set contrac 0.07
set dialt1 0.4
set dialt2 2
set liquefac1 10
set liquefac2 0.01
set liquefac3 1
set e 0.7
# Number of yield surfaces
set noYieldSurf 20
# General Parameters
set g 9.81

puts "End Of Define Parameters"

#-------------------------------------------------------------------------------------------
# 2. DEFINE NODES FOR SOIL ELEMENTS
#-------------------------------------------------------------------------------------------

model BasicBuilder -ndm 2 -ndf 2 ;

file mkdir Input_Files
set NodeInfo1 [open Input_Files/nodeinfo.dat w]

for {set j 1} {$j <=$NNodeY1} {incr j 1} {
for {set i 1} {$i <= $NNodeX} {incr i 1} {
set NodeN1 [expr (($j-1)*$NNodeX)+$i]
set XCoor($NodeN1) [expr ($i-1)*$DeltaX]
set YCoor($NodeN1) [expr ($j-1)*$DeltaY]
node $NodeN1 $XCoor($NodeN1) $YCoor($NodeN1)
puts $NodeInfo1 "node $NodeN1 $XCoor($NodeN1) $YCoor($NodeN1)"

}
};
close $NodeInfo1
puts "Finished creating all soil nodes..."

#-------------------------------------------------------------------------------------------
# 3. DEFINE BOUNDARY CONDITIONS
#-------------------------------------------------------------------------------------------

set BoundryCondition1 [open Input_Files/BoundryCondition.dat w]

fix 1 1 1
fix $NNodeX 1 1
puts "$NNodeX"

for {set i 2} {$i<=$NNodeX-1} {incr i 1} {
fix $i 0 1
puts $BoundryCondition1 "node fix $i 0 1"
};

for {set j 2} {$j <= $NNodeY1} {incr j 1} {
equalDOF [expr (($j-1)*$NNodeX)+1] [expr (($j-1)*$NNodeX)+$NNodeX] 1
puts $BoundryCondition1 "equal Dof [expr (($j-1)*$NNodeX)+1] [expr (($j-1)*$NNodeX)+$NNodeX] 1"
};

close $BoundryCondition1
puts "End Of Define Boundry Conditions"

#-------------------------------------------------------------------------------------------
# 4. DEFINE SOIL MATERIALS
#-------------------------------------------------------------------------------------------
# Define Soil Materials

nDMaterial PressureDependMultiYield 1 2 $rho $refShearModul $refBulkModul $frictionAng $peakShearStra $refPress $pressDependCoe $PTAng $contrac $dialt1 $dialt2 $liquefac1 $liquefac2 $liquefac3 $noYieldSurf

puts "Finished creating all soil materials..."

#-------------------------------------------------------------------------------------------
# 5. DEFINE SOIL ELEMENTS
#-----------------------------------------------------------------------------------------
set EleInfo1 [open Input_Files/eleinfo.dat w]

set wgtX 0.0
set wgtY [expr -9.81*$rho]

for {set j 1} {$j <= $NeleY1} {incr j 1} {
for {set i 1} {$i <= $NeleX} {incr i 1} {
set ElemN1 [expr (($j-1)*$NeleX)+$i]
set Node1 [expr ($j-1)*$NNodeX+$i]
set Node2 [expr $Node1+1]
set Node3 [expr $j*$NNodeX+($i+1)]
set Node4 [expr $Node3-1]
element quad $ElemN1 $Node1 $Node2 $Node3 $Node4 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY
puts $EleInfo1 "Element Quad $ElemN1 $Node1 $Node2 $Node3 $Node4 thick=1.0 rho=$rho unitWeightX=$wgtX unitWeightY=$wgtY"
}
};

close $EleInfo1
puts "Finished creating all soil elements..."

file mkdir Output_files

updateMaterialStage -material 1 -stage 0

recorder Element -file Output_files/Gstress.txt -time -eleRange 1 4 material 1 stress

#-------------------------------------------------------------------------------------------
# 6. GRAVITY APPLICATION
#-----------------------------------------------------------------------------------------

numberer $numberer
eval "system $solver"
test NormDispIncr 1.00e-002 30 0
algorithm ModifiedNewton
constraints Penalty 1.e18 1.e18
#integrator LoadControl 1 1 1 1
#analysis Static
integrator Newmark 1.5 1.
analysis Transient

set ok [analyze 5 5.00e+004]
if {$ok != 0} { return $ok }
puts "First run done."

# switch material stage from elastic (gravity) to plastic ------- 2nd run

updateMaterialStage -material 1 -stage 1

# plastic gravity loading
set ok [analyze 5 5.00e+004]
if {$ok != 0} { return $ok }
puts "Finished with plastic gravity analysis..."

#-------------------------------------------------------------------------------------------
# 7. CREATE NODES AND ELEMENTS FOR TOP SOIL
#-------------------------------------------------------------------------------------------

#----NODES
set NodeInfo2 [open Input_Files/nodeinfo2.dat w]
for {set j 1} {$j <=$NNodeY2} {incr j 1} {
for {set i 1} {$i <= $NNodeX} {incr i 1} {
set NodeN2 [expr $NodeN1+($j-1)*$NNodeX+$i]
set XCoor($NodeN2) [expr ($i-1)*$DeltaX]
set YCoor($NodeN2) [expr ($j-1)*$DeltaY+($Width+($DeltaY*$NeleY1))]
node $NodeN2 $XCoor($NodeN2) $YCoor($NodeN2)
puts $NodeInfo2 "node $NodeN2 $XCoor($NodeN2) $YCoor($NodeN2)"

}
};
close $NodeInfo2

#----BOUNDRY CONDITION
set BoundryCondition2 [open Input_Files/BoundryCondition2.dat w]
for {set j 1} {$j <= $NNodeY2} {incr j 1} {
equalDOF [expr ($NNodeX*$NNodeY1)+(($j-1)*$NNodeX)+1] [expr ($NNodeX*$NNodeY1)+(($j-1)*$NNodeX)+$NNodeX] 1 2
puts $BoundryCondition2 "equal Dof [expr ($NNodeX*$NNodeY1)+(($j-1)*$NNodeX)+1] [expr ($NNodeX*$NNodeY1)+(($j-1)*$NNodeX)+$NNodeX] 1 2"
};
close $BoundryCondition2

#----TOP SOIL ELEMENTS
set EleInfo2 [open Input_Files/EleInfo2.dat w]
for {set j 1} {$j <= $NeleY2} {incr j 1} {
for {set i 1} {$i <= $NeleX} {incr i 1} {
set ElemN2 [expr ($NeleX*$NeleY1)+(($j-1)*$NeleX)+$i]
set Node1 [expr ($NNodeX*$NNodeY1)+(($j-1)*$NNodeX)+$i]
set Node2 [expr $Node1+1]
set Node3 [expr ($NNodeX*$NNodeY1)+(($j-1)*$NNodeX)+$NNodeX+1+$i]
set Node4 [expr $Node3-1]
element quad $ElemN2 $Node1 $Node2 $Node3 $Node4 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY
puts $EleInfo2 "Element Quad $ElemN2 $Node1 $Node2 $Node3 $Node4 thick=5.0 rho=$rho unitWeightX=$wgtX unitWeightY=$wgtY"
}
};

close $EleInfo2
puts "Finished creating Nodes and Elements For Top Soil... "

#-------------------------------------------------------------------------------------------
# 8. CREATE LAGRANGE MULTIPLIER NODES FOR BEAM CONTACT ELEMENTS
#-------------------------------------------------------------------------------------------

set LagrangeNodes [open Input_Files/LagrangeNodes.dat w]
for {set k 1} {$k <= [expr 2*$NNodeX]} {incr k 1} {
set LNodes [expr $NodeN2+$k]
set X($LNodes) 0.0
set Y($LNodes) [expr ($NeleY1*$DeltaY)+($Width/2)]
node $LNodes $X($LNodes) $Y($LNodes)
puts $LagrangeNodes "$LNodes $X($LNodes) $Y($LNodes)"

};

close $LagrangeNodes
puts "Finished creating all -ndf 2 nodes"

#-------------------------------------------------------------------------------------------
# 9. DEFINE BEAM NODES
#-----------------------------------------------------------------------------------------

model BasicBuilder -ndm 2 -ndf 3 ;

set NodeInfoP [open Input_Files/nodeinfoP.dat w]

for {set k 1} {$k<=$NNodeX+1} {incr k 1} {
set ijNode [expr $LNodes+$k]
set X($ijNode) [expr (($k-1)*$DeltaX)-($DeltaX/2)]
set Y($ijNode) [expr ($NeleY1*$DeltaY)+($Width/2)]
node $ijNode $X($ijNode) $Y($ijNode)
puts $NodeInfoP "Pipeline Nodes $ijNode $X($ijNode) $Y($ijNode)"

}

close $NodeInfoP
puts "Finished Creating Beam Nodes..."

#-------------------------------------------------------------------------------------------
# 10. DEFINE BEAM MATERIAL
#----------------------------------------------------------------------------------------

#----PIPE MATERIAL

set Fy 358000 ;# kPa
set Eb 2.0e8 ;# kPa
set b 0.01 ;# strain hardening ratio
set R0 18
set cR1 0.925
set cR2 0.15

uniaxialMaterial Steel02 4 $Fy $Eb $b $R0 $cR1 $cR2

#----PIPE SECTION

section Fiber 1 {

set numSubdivCirc 12
set numSubdivRad 4
set yCenter 0.0
set zCenter 0.0
set intRad 0.09685 ;# internal radious
set extRad 0.10955 ;# external radious
set startAng 0.0
set endAng [expr 2*$pi]
patch circ 4 $numSubdivCirc $numSubdivRad $yCenter $zCenter $intRad $extRad $startAng $endAng

}

#----GEOMETRIC TRANSFORMATION

geomTransf Linear 1

puts "Finished Creating Beam Material..."

#-------------------------------------------------------------------------------------------
# 11. DEFINE BEAM ELEMENTS
#----------------------------------------------------------------------------------------

set EleInfoP [open Input_Files/eleinfoP.dat w]

for {set p 1} {$p<=$NNodeX} {incr p 1} {
set EleP [expr $ElemN2+$p]
set numIntgrPts 5
set massDens 64.58 ;#kg/m
element nonlinearBeamColumn $EleP [expr $LNodes+$p] [expr $LNodes+$p+1] $numIntgrPts 1 1 $massDens

puts $EleInfoP "Nonlinear Beam Elements $EleP iNode=[expr $LNodes+$p] jNode=[expr $LNodes+$p+1]"
}

close $EleInfoP
puts "Finished Define Beam Elements..."

#-------------------------------------------------------------------------------------------
# 12. DEFINE CONTACT MATERIAL FOR BEAM CONTACT ELEMENTS
#----------------------------------------------------------------------------------------

# two-dimensional contact material

set f 0.8 ;#Firiction factor for Rough Steel
set D [expr $extRad*2] ;#Outside Diameter of Pipe
set c 0.0 ;#cohesive soil
set Hs [expr $V2+$extRad] ;#Depth of soil above the center of the pipe
set K0 [expr 1-sin($phi)]
set deltaprim [expr $f*$frictionAng]
set deltaprim2 [expr ($deltaprim*$pi)/180] ;#Radian
set mu [expr tan($deltaprim2)] ;#nterface frictional coefficient
set Gs 2.66 ;#Specific Gravity Of Soil
set GammaW 9.81 ;#Kg/m^3
set EffectiveGamma [expr ($Gs-1)*$GammaW/(1+$e)]
set tu [expr $pi*$D*$Hs*$EffectiveGamma*((1+$K0)/2)*tan($deltaprim2)]
set Deltamm 0.003 ;#m For Dense Sand
set Kx [expr $tu/$Deltamm] ;#interface stiffness parameter
set t 0.0
puts "Kx=$Kx" ;#interface tensile strength

#Contact material 6
set Nq 32
set Ngamma 22
set Qd [expr ($Nq*$EffectiveGamma*$Hs*$D)+($Ngamma*$rho*$g*$D*$D/2)]
set Deltaqd [expr 0.1*$D]
set Kz [expr $Qd/$Deltaqd]
puts "Kz=$Kz"

nDMaterial ContactMaterial2D 5 $mu $Kz 0.0 0.0

puts "Finished creating all contact materials..."

#-----------------------------------------------------------------------------------------
# 13. CREATE BEAM CONTACT ELEMENTS
#-----------------------------------------------------------------------------------------

set ContactBeam [open Input_Files/ContactBeam.dat w]

# set gap and force tolerances for beam contact elements

set gTol 1.0e-10
set fTol 1.0e-10
set cFlag 0

# define beam contact elements

for {set M 1} {$M <= $NNodeX} {incr M 1} {

set EleBC1 [expr $EleP+$M]
set Snode1 [expr ($NNodeY1-1)*$NNodeX+$M]
set LagrangeNode1 [expr $NodeN2+$M]
element BeamContact2D $EleBC1 [expr $LNodes+$M] [expr $LNodes+$M+1] $Snode1 $LagrangeNode1 5 $Width $gTol $fTol
puts $ContactBeam "EleBC1=$EleBC1 iNode=[expr $LNodes+$M] jNode=[expr $LNodes+$M+1] Snode1=$Snode1 LagrangeNode1=$LagrangeNode1 "
}

for {set E 1} {$E <= $NNodeX} {incr E 1} {

set EleBC2 [expr $EleBC1+$E]
set Snode2 [expr ($NNodeX*$NNodeY1)+$E]
set LagrangeNode2 [expr $LagrangeNode1+$E]
element BeamContact2D $EleBC2 [expr $LNodes+$E] [expr $LNodes+$E+1] $Snode2 $LagrangeNode2 5 $Width $gTol $fTol
puts $ContactBeam "EleBC2=$EleBC2 iNode=[expr $LNodes+$E] jNode=[expr $LNodes+$E+1] Snode2=$Snode2 LagrangeNode2=$LagrangeNode2"
}

setParameter -value 1 -eleRange [expr $EleP+1] [expr $EleP+(2*$NNodeX)] friction

close $ContactBeam
puts "Finished creating all contact Beam..."

#-------------------------------------------------------------------------------------------
# 14. CREATE RECORDERS
#------------------------------------------------------------------------------------------
file mkdir Output_files

for {set R 1} {$R <= $NNodeX} {incr R 1} {
set ele [expr $ElemN2+$R]
set sec $R
set y [expr ($R-1)*$DeltaX]
set z [expr ($NeleY1*$DeltaY)+($Width/2)]
recorder Element -file Output_files/stressStrain.out –time -ele $ele section $sec fiber $y $z stressStrain
};

# View Model
recorder display "Model" 10 10 1300 700 -wipe
prp 0 0 50
vup 0 1 0
vpn 0 0 1
display 1 2 10

puts "Finished creating Time History recorders..."

#-------------------------------------------------------------------------------------------
# 15. DYNAMIC ANALYSIS
#-------------------------------------------------------------------------------------------

model BasicBuilder -ndm 2 -ndf 2 ;

set dt 0.01
set Tmax 40.95

set Gaccel "Series -dt $dt -filePath 7_1.txt -factor 9.81" ;
pattern UniformExcitation 4 1 -accel "Series -dt $dt -filePath 7_1.txt -factor 9.81 -startTime [getTime]";

puts "Dynamic loading created..."

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

#---ANALYSIS PARAMETERS
# Newmark parameters
set gamma 0.5
set beta 0.25

rayleigh $a0 $a1 0.0 0.0

constraints Transformation ;
numberer RCM ;
system SparseGeneral ;
test EnergyIncr 1e-1 200 1 ;
algorithm KrylovNewton ;
integrator TRBDF2 ;
analysis VariableTransient ;
analyze [expr int($duration/$dt)] $dt ;

puts "Finished with dynamic analysis...";