Handling "Transformation" Constraints

A forum dedicated to users with questions regarding soil materials and elements.

forum currently locked

Moderator: Moderators

Locked
angelo.lambrughi
Posts: 4
Joined: Fri Mar 04, 2011 3:48 am
Location: SOIL srl

Handling "Transformation" Constraints

Post by angelo.lambrughi »

Dear all,

I've been working in the past few days on a 2D numerical model of a submerged slope.

I began from the very useful example problem available within the OPENSEES documentation ("Dynamic Effective Stress Analysis of a Slope"). However, in that problem fixed null values are assigned to pore pressures along the "mudline" of the slope.

In order to apply proper fixed hydrostatic values of pore pressure along the boundary, I then followed the procedure suggested within the other example "Simulating a Centrifuge Test", which requires the enforcement of "Transformation" constraints (instead of the "Penalty" method, adopted in the slope example problem).

Unfortunately, I got a kind of conflict between the "equaldof" conditions, which are assigned to obtain "periodic boundary conditions" and the "Transformation" constraints. I'm aware of the restriction one should comply with when "Transformation" constraints are used (a "master" node in a "equaldof" condition cannot play the role pf "slave" node in another condition, for instance).

In any case, I wrote basic models in order to better investigate my problem.

While a very basic model made by a single "9_4_quadUP" element seems works properly (though at a very slow pace), I was not able to bring the end the analysis of an eight "9_4_quadUP" system.

In both models, a fixed value of 30 kPa is assigned at the top face for pore pressure. For total equilibrium, equivalent (and consistent) nodal loads are applied at the same surface.

For the single element model, proper values of pore pressure and effective stresses are obtained (though for this very simple case no periodic boundary conditions can be enforced). For the eight element models, instead, I cannot even reach convergence during the gravity stage.

I'm almost sure that my problems are due to improper assignment of "equaldof" conditions.

Anyone has any suggestions?

In the following the two basic models codes are posted.

----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
MODEL 1 - SINGLE ELEMENT
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

#-----------------------------------------------------------------------------------------
# 1. CREATE PORE PRESSURE NODES AND FIXITIES
#-----------------------------------------------------------------------------------------

model BasicBuilder -ndm 2 -ndf 3

# define pore pressure nodes:

node 1 0.000 0
node 2 2.000 0.000
node 3 2.000 2.000
node 4 0.000 2.000

set mNodeInfo [open nodeInfo.dat w]

puts $mNodeInfo " 1 0.000 0.000"
puts $mNodeInfo " 2 2.000 0.000"
puts $mNodeInfo " 3 2.000 2.000"
puts $mNodeInfo " 4 0.000 2.000"

close $mNodeInfo

puts "Finished creating all -ndf 3 nodes..."

# define fixities for pore pressure nodes:

fix 1 0 1 0
fix 2 0 1 0
fix 4 0 0 1
fix 3 0 0 1

set Timeseries03 "Constant -factor 30"

set Timeseries04 "Constant -factor 30"

set Timeseries00 "Constant -factor 0"

pattern MultipleSupport 3 {

groundMotion 100003 Plain -vel $Timeseries03

imposedMotion 3 3 100003

groundMotion 100004 Plain -vel $Timeseries04

imposedMotion 4 3 100004

}

puts "Finished creating all -ndf 3 boundary conditions..."

# define equal degrees of freedom for pore pressure nodes

equalDOF 1 2 1
equalDOF 4 3 1

puts "Finished creating equalDOF for pore pressure nodes..."

#-----------------------------------------------------------------------------------------
# 2. CREATE INTERIOR NODES AND FIXITIES
#-----------------------------------------------------------------------------------------
model BasicBuilder -ndm 2 -ndf 2

# define interior nodes

node 5 1.000 0.000
node 6 2.000 1.00
node 7 1.000 2.00
node 8 0.000 1.00
node 9 1.000 1.0

puts "Finished creating all -ndf 2 nodes..."


# define fixities for interior nodes:

fix 5 0 1

puts "Finished creating all -ndf 2 boundary conditions..."

# define equal degrees of freedom which have not yet been defined

#equalDOF 1 5 1
#equalDOF 8 9 1 2
equalDOF 8 6 1 2
#equalDOF 4 7 1 2

puts "Finished creating equalDOF constraints..."


puts "Finished creating equalDOF for base..."

#-----------------------------------------------------------------------------------------
# 3. CREATE SOIL MATERIALS
#-----------------------------------------------------------------------------------------

nDMaterial PressureDependMultiYield02 1 2 1.8 14580 24300 32 3 \
14 0.5 30 0.087 0.18 0 \
0 5 5.0 3.0 1.0 \
0.0 0.65 0.9 0.02 0.7 101.0

set thick1 1.0
set xWgt1 0.00
set yWgt1 -9.81
set uBulk1 5.3e5
set hPerm1 1.0e-4
set vPerm1 1.0e-4

#-----------------------------------------------------------------------------------------
# 4. CREATE SOIL ELEMENTS
#-----------------------------------------------------------------------------------------

# permeabilities are initial set at 1.0 m/s for gravity analysis, values are updated after gravity

element 9_4_QuadUP 1 1 2 3 4 5 6 7 8 9 $thick1 1 $uBulk1 1.0 1.0 1.0 $xWgt1 $yWgt1

puts "Finished creating all soil elements..."

#-----------------------------------------------------------------------------------------
# 5. LYSMER DASHPOT
#-----------------------------------------------------------------------------------------

# define dashpot nodes
node 16 0.000 0.000
node 17 0.000 0.000

# define fixities for dashpot nodes
fix 16 1 1
fix 17 0 1

# define equal DOF for dashpot and base soil node
equalDOF 1 17 1
puts "Finished creating dashpot nodes and boundary conditions..."

# define dashpot material
set baseArea 1.00
set dashpotCoeff 664.20
uniaxialMaterial Viscous 2 [expr $dashpotCoeff*$baseArea] 1

# define dashpot element
element zeroLength 3 16 17 -mat 2 -dir 1
puts "Finished creating dashpot material and element..."

#-----------------------------------------------------------------------------------------
# 6. DEFINE NODAL MASSES FOR MODELING WATER
#-----------------------------------------------------------------------------------------

# define nodal masses for 3 dof nodes
model BasicBuilder -ndm 2 -ndf 3

mass 3 0 1.02 0

mass 4 0 1.02 0

puts "Finished creating -ndf 3 nodal masses..."

# define nodal masses for 2 dof nodes
model BasicBuilder -ndm 2 -ndf 2

mass 7 0 4.08

puts "Finished creating -ndf 2 nodal masses..."

#-----------------------------------------------------------------------------------------
# 7. CREATE GRAVITY RECORDERS
#-----------------------------------------------------------------------------------------

# create list for pore pressure nodes
set nodeList3 {}
set channel [open "nodeInfo.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

recorder Element -file Gstress1.out -time -eleRange 1 1 material 1 stress
recorder Element -file Gstress2.out -time -eleRange 1 1 material 2 stress
recorder Element -file Gstress3.out -time -eleRange 1 1 material 3 stress
recorder Element -file Gstress4.out -time -eleRange 1 1 material 4 stress
recorder Element -file Gstress5.out -time -eleRange 1 1 material 5 stress
recorder Element -file Gstress6.out -time -eleRange 1 1 material 6 stress
recorder Element -file Gstress7.out -time -eleRange 1 1 material 7 stress
recorder Element -file Gstress8.out -time -eleRange 1 1 material 8 stress
recorder Element -file Gstress9.out -time -eleRange 1 1 material 9 stress
recorder Element -file Gstrain1.out -time -eleRange 1 1 material 1 strain
recorder Element -file Gstrain2.out -time -eleRange 1 1 material 2 strain
recorder Element -file Gstrain3.out -time -eleRange 1 1 material 3 strain
recorder Element -file Gstrain4.out -time -eleRange 1 1 material 4 strain
recorder Element -file Gstrain5.out -time -eleRange 1 1 material 5 strain
recorder Element -file Gstrain6.out -time -eleRange 1 1 material 6 strain
recorder Element -file Gstrain7.out -time -eleRange 1 1 material 7 strain
recorder Element -file Gstrain8.out -time -eleRange 1 1 material 8 strain
recorder Element -file Gstrain9.out -time -eleRange 1 1 material 9 strain

recorder Element -file G_Gauss1.out -time -eleRange 1 1 material 1 gauss

puts "Finished creating gravity recorders..."

#-----------------------------------------------------------------------------------------
# 8. CREATE FILES FOR POSTPROCESSING IN GiD
#-----------------------------------------------------------------------------------------

set meshFile [open renameMe.flavia.msh w]
puts $meshFile "MESH 94quad dimension 2 ElemType Quadrilateral Nnode 4"
puts $meshFile "Coordinates"
puts $meshFile "#node_number coord_x coord_y"
puts $meshFile " 1 0.000 0.000"
puts $meshFile " 2 2.000 0.000"
puts $meshFile " 3 2.000 2.000"
puts $meshFile " 4 0.000 2.000"

puts $meshFile "end coordinates"

puts $meshFile "Elements"
puts $meshFile "# element node1 node2 node3 node4"
puts $meshFile " 1 1 2 3 4"

puts $meshFile "end elements"

close $meshFile

set eleFile [open elementInfo.dat w]

puts $eleFile " 1 1 2 3 4"

close $eleFile

#-----------------------------------------------------------------------------------------
# 9. DEFINE ANALYSIS PARAMETERS
#-----------------------------------------------------------------------------------------

#---GROUND MOTION PARAMETERS
# time step in ground motion record

set motionDT 0.01

# number of steps in ground motion record

set motionSteps 2987

#---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 369.0
# element size (m)

set eleSize 2
# duration of ground motion (s)
set duration [expr $motionDT*$motionSteps]
# trial analysis time step
set kTrial [expr $eleSize/(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
#-----------------------------------------------------------------------------------------

model BasicBuilder -ndm 2 -ndf 3

# loading object
pattern Plain 20 Constant {

load 3 0 -10 0

load 4 0 -10 0

}

model BasicBuilder -ndm 2 -ndf 2

pattern Plain 30 Constant {

load 7 0 -40

}

model BasicBuilder -ndm 2 -ndf 3

# update materials to ensure elastic behavior

updateMaterialStage -material 1 -stage 0

set gamma 1.6

# create the SOE, ConstraintHandler, Integrator, Algorithm and Numberer

integrator Newmark $gamma [expr pow($gamma+0.5, 2)/4] 0.00 0.0 0.00 0.0

test EnergyIncr 1.0e-8 400 1;

constraints Transformation

algorithm Newton

numberer RCM

system ProfileSPD

analysis Transient

set startT [clock seconds]
analyze 10 5.0e2
analyze 10 5.0e3
puts "Finished with elastic gravity analysis..."

# update materials to consider plastic behavior

updateMaterialStage -material 1 -stage 1

# plastic gravity loading
analyze 10 5.0e-3
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<=1} {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
updateParameter 10001 $vPerm1
updateParameter 10002 $hPerm1

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

recorder Element -file stress1.out -time -dT $recDT -eleRange 1 1 material 1 stress
recorder Element -file stress2.out -time -dT $recDT -eleRange 1 1 material 2 stress
recorder Element -file stress3.out -time -dT $recDT -eleRange 1 1 material 3 stress
recorder Element -file stress4.out -time -dT $recDT -eleRange 1 1 material 4 stress
recorder Element -file stress5.out -time -dT $recDT -eleRange 1 1 material 5 stress
recorder Element -file stress6.out -time -dT $recDT -eleRange 1 1 material 6 stress
recorder Element -file stress7.out -time -dT $recDT -eleRange 1 1 material 7 stress
recorder Element -file stress8.out -time -dT $recDT -eleRange 1 1 material 8 stress
recorder Element -file stress9.out -time -dT $recDT -eleRange 1 1 material 9 stress
recorder Element -file strain1.out -time -dT $recDT -eleRange 1 1 material 1 strain
recorder Element -file strain2.out -time -dT $recDT -eleRange 1 1 material 2 strain
recorder Element -file strain3.out -time -dT $recDT -eleRange 1 1 material 3 strain
recorder Element -file strain4.out -time -dT $recDT -eleRange 1 1 material 4 strain
recorder Element -file strain5.out -time -dT $recDT -eleRange 1 1 material 5 strain
recorder Element -file strain6.out -time -dT $recDT -eleRange 1 1 material 6 strain
recorder Element -file strain7.out -time -dT $recDT -eleRange 1 1 material 7 strain
recorder Element -file strain8.out -time -dT $recDT -eleRange 1 1 material 8 strain
recorder Element -file strain9.out -time -dT $recDT -eleRange 1 1 material 9 strain
puts "Finished creating all recorders..."

#-----------------------------------------------------------------------------------------
# 13. DYNAMIC ANALYSIS
#-----------------------------------------------------------------------------------------

model BasicBuilder -ndm 2 -ndf 3

# define constant scaling factor for applied velocity
set cFactor [expr $baseArea*$dashpotCoeff]

# define velocity time history file
set velocityFile velocityHistory.out

# timeseries object for force history
set mSeries "Path -dt $motionDT -filePath $velocityFile -factor $cFactor"

# loading object
pattern Plain 10 $mSeries {
load 1 1.0 0.0 0.0
}
puts "Dynamic loading created..."

set gamma 1.6

# create the SOE, ConstraintHandler, Integrator, Algorithm and Numberer

integrator Newmark $gamma [expr pow($gamma+0.5, 2)/4] 0.00 0.0 0.00 0.0

#test EnergyIncr 1.0e-8 400 1;

test NormDispIncr 1.0e-3 400 1

constraints Transformation

algorithm Newton

numberer RCM

system ProfileSPD

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]
puts "curTime: $curTime"
set curStep [expr $curTime/$dT]
puts "curStep: $curStep"
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/$dT]
puts "curStep: $curStep"
set remStep [expr int(($remStep-$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

----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
MODEL 2 - EIGHT ELEMENTS
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

#-----------------------------------------------------------------------------------------
# 1. CREATE PORE PRESSURE NODES AND FIXITIES
#-----------------------------------------------------------------------------------------

model BasicBuilder -ndm 2 -ndf 3

# define pore pressure nodes:

node 1 0.000 0
node 2 2.000 0.000
node 3 2.000 2.000
node 4 0.000 2.000
node 10 2 4
node 11 0 4
node 16 4 0
node 17 4 2
node 22 4 4
node 26 6 0
node 27 6 2
node 32 6 4
node 36 8 0
node 37 8 2
node 42 8 4

set mNodeInfo [open nodeInfo.dat w]

puts $mNodeInfo " 1 0.000 0.000"
puts $mNodeInfo " 2 2.000 0.000"
puts $mNodeInfo " 3 2.000 2.000"
puts $mNodeInfo " 4 0.000 2.000"
puts $mNodeInfo " 10 2 4"
puts $mNodeInfo " 11 0 4"
puts $mNodeInfo " 16 4 0"
puts $mNodeInfo " 17 4 2"
puts $mNodeInfo " 22 4 4"
puts $mNodeInfo " 26 6 0"
puts $mNodeInfo " 27 6 2"
puts $mNodeInfo " 32 6 4"
puts $mNodeInfo " 36 8 0"
puts $mNodeInfo " 37 8 2"
puts $mNodeInfo " 42 8 4"

close $mNodeInfo

puts "Finished creating all -ndf 3 nodes..."

# define fixities for pore pressure nodes:

fix 1 0 1 0
fix 2 0 1 0
fix 16 0 1 0
fix 26 0 1 0
fix 36 0 1 0

fix 11 0 0 1
fix 10 0 0 1
fix 22 0 0 1
fix 32 0 0 1
fix 42 0 0 1

set Timeseries11 "Constant -factor 30"

set Timeseries10 "Constant -factor 30"

set Timeseries22 "Constant -factor 30"

set Timeseries32 "Constant -factor 30"

set Timeseries42 "Constant -factor 30"

set Timeseries00 "Constant -factor 0"

pattern MultipleSupport 3 {

groundMotion 100011 Plain -vel $Timeseries11

imposedMotion 11 3 100011

groundMotion 100010 Plain -vel $Timeseries10

imposedMotion 10 3 100010

groundMotion 100022 Plain -vel $Timeseries22

imposedMotion 22 3 100022

groundMotion 100032 Plain -vel $Timeseries32

imposedMotion 32 3 100032

groundMotion 100042 Plain -vel $Timeseries42

imposedMotion 42 3 100042

}

puts "Finished creating all -ndf 3 boundary conditions..."

# define equal degrees of freedom for pore pressure nodes

equalDOF 1 2 1

equalDOF 4 3 1

equalDOF 11 10 1

equalDOF 27 37 1

equalDOF 32 42 1

puts "Finished creating equalDOF for pore pressure nodes..."

#-----------------------------------------------------------------------------------------
# 2. CREATE INTERIOR NODES AND FIXITIES
#-----------------------------------------------------------------------------------------
model BasicBuilder -ndm 2 -ndf 2

# define interior nodes

node 5 1 0
node 6 2 1
node 7 1 2
node 8 0 1
node 9 1 1
node 12 2 3
node 13 1 4
node 14 0 3
node 15 1 3
node 18 3 0
node 19 4 1
node 20 3 2
node 21 3 1
node 23 4 3
node 24 3 4
node 25 3 3
node 28 5 0
node 29 6 1
node 30 5 2
node 31 5 1
node 33 6 3
node 34 5 4
node 35 5 3
node 38 7 0
node 39 8 1
node 40 7 2
node 41 7 1
node 43 8 3
node 44 7 4
node 45 7 3

puts "Finished creating all -ndf 2 nodes..."


# define fixities for interior nodes:

fix 5 0 1

fix 18 0 1

fix 28 0 1

fix 38 0 1

puts "Finished creating all -ndf 2 boundary conditions..."

# define equal degrees of freedom which have not yet been defined

equalDOF 8 6 1

equalDOF 14 12 1

equalDOF 29 39 1

equalDOF 33 43 1


puts "Finished creating equalDOF constraints..."

equalDOF 1 5 1

equalDOF 1 18 1

equalDOF 1 16 1

equalDOF 1 28 1

equalDOF 1 26 1

equalDOF 1 38 1

equalDOF 1 36 1

puts "Finished creating equalDOF for base..."

#-----------------------------------------------------------------------------------------
# 3. CREATE SOIL MATERIALS
#-----------------------------------------------------------------------------------------

nDMaterial PressureDependMultiYield02 1 2 1.8 14580 24300 32 3 \
14 0.5 30 0.087 0.18 0 \
0 5 5.0 3.0 1.0 \
0.0 0.65 0.9 0.02 0.7 101.0

set thick1 1.0
set xWgt1 0.00
set yWgt1 -9.81
set uBulk1 5.3e5
set hPerm1 1.0e-4
set vPerm1 1.0e-4

nDMaterial PressureDependMultiYield02 2 2 1.8 14580 24300 32 3 \
14 0.5 30 0.087 0.18 0 \
0 5 5.0 3.0 1.0 \
0.0 0.65 0.9 0.02 0.7 101.0

set thick2 10000
set xWgt2 0.00
set yWgt2 -9.81
set uBulk2 5.3e5
set hPerm2 1.0e-4
set vPerm2 1.0e-4

#-----------------------------------------------------------------------------------------
# 4. CREATE SOIL ELEMENTS
#-----------------------------------------------------------------------------------------

# permeabilities are initial set at 1.0 m/s for gravity analysis, values are updated after gravity

element 9_4_QuadUP 1 1 2 3 4 5 6 7 8 9 $thick2 2 $uBulk2 1.0 1.0 1.0 $xWgt2 $yWgt2

element 9_4_QuadUP 2 4 3 10 11 7 12 13 14 15 $thick2 2 $uBulk2 1.0 1.0 1.0 $xWgt2 $yWgt2

element 9_4_QuadUP 3 2 16 17 3 18 19 20 6 21 $thick1 1 $uBulk1 1.0 1.0 1.0 $xWgt1 $yWgt1

element 9_4_QuadUP 4 3 17 22 10 20 23 24 12 25 $thick1 1 $uBulk1 1.0 1.0 1.0 $xWgt1 $yWgt1

element 9_4_QuadUP 5 16 26 27 17 28 29 30 19 31 $thick1 1 $uBulk1 1.0 1.0 1.0 $xWgt1 $yWgt1

element 9_4_QuadUP 6 17 27 32 22 30 33 34 23 35 $thick1 1 $uBulk1 1.0 1.0 1.0 $xWgt1 $yWgt1

element 9_4_QuadUP 7 26 36 37 27 38 39 40 29 41 $thick2 2 $uBulk2 1.0 1.0 1.0 $xWgt2 $yWgt2

element 9_4_QuadUP 8 27 37 42 32 40 43 44 33 45 $thick2 2 $uBulk2 1.0 1.0 1.0 $xWgt2 $yWgt2

puts "Finished creating all soil elements..."

#-----------------------------------------------------------------------------------------
# 5. LYSMER DASHPOT
#-----------------------------------------------------------------------------------------

# define dashpot nodes
node 46 0.000 0.000
node 47 0.000 0.000

# define fixities for dashpot nodes
fix 46 1 1
fix 47 0 1

# define equal DOF for dashpot and base soil node
equalDOF 1 47 1
puts "Finished creating dashpot nodes and boundary conditions..."

# define dashpot material
set baseArea 40004
set dashpotCoeff 664.20
uniaxialMaterial Viscous 3 [expr $dashpotCoeff*$baseArea] 1

# define dashpot element
element zeroLength 9 46 47 -mat 3 -dir 1
puts "Finished creating dashpot material and element..."

#-----------------------------------------------------------------------------------------
# 6. DEFINE NODAL MASSES FOR MODELING WATER
#-----------------------------------------------------------------------------------------

# define nodal masses for 3 dof nodes
model BasicBuilder -ndm 2 -ndf 3

mass 11 0 1.02 0

mass 10 0 2.04 0

mass 22 0 2.04 0

mass 32 0 2.04 0

mass 42 0 1.02 0

puts "Finished creating -ndf 3 nodal masses..."

# define nodal masses for 2 dof nodes
model BasicBuilder -ndm 2 -ndf 2

mass 13 0 4.08

mass 24 0 4.08

mass 34 0 4.08

mass 44 0 4.08

puts "Finished creating -ndf 2 nodal masses..."

#-----------------------------------------------------------------------------------------
# 7. CREATE GRAVITY RECORDERS
#-----------------------------------------------------------------------------------------

# create list for pore pressure nodes
set nodeList3 {}
set channel [open "nodeInfo.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

recorder Element -file Gstress1.out -time -eleRange 1 8 material 1 stress
recorder Element -file Gstress2.out -time -eleRange 1 8 material 2 stress
recorder Element -file Gstress3.out -time -eleRange 1 8 material 3 stress
recorder Element -file Gstress4.out -time -eleRange 1 8 material 4 stress
recorder Element -file Gstress5.out -time -eleRange 1 8 material 5 stress
recorder Element -file Gstress6.out -time -eleRange 1 8 material 6 stress
recorder Element -file Gstress7.out -time -eleRange 1 8 material 7 stress
recorder Element -file Gstress8.out -time -eleRange 1 8 material 8 stress
recorder Element -file Gstress9.out -time -eleRange 1 8 material 9 stress
recorder Element -file Gstrain1.out -time -eleRange 1 8 material 1 strain
recorder Element -file Gstrain2.out -time -eleRange 1 8 material 2 strain
recorder Element -file Gstrain3.out -time -eleRange 1 8 material 3 strain
recorder Element -file Gstrain4.out -time -eleRange 1 8 material 4 strain
recorder Element -file Gstrain5.out -time -eleRange 1 8 material 5 strain
recorder Element -file Gstrain6.out -time -eleRange 1 8 material 6 strain
recorder Element -file Gstrain7.out -time -eleRange 1 8 material 7 strain
recorder Element -file Gstrain8.out -time -eleRange 1 8 material 8 strain
recorder Element -file Gstrain9.out -time -eleRange 1 8 material 9 strain

puts "Finished creating gravity recorders..."

#-----------------------------------------------------------------------------------------
# 8. CREATE FILES FOR POSTPROCESSING IN GiD
#-----------------------------------------------------------------------------------------

set meshFile [open renameMe.flavia.msh w]
puts $meshFile "MESH 94quad dimension 2 ElemType Quadrilateral Nnode 4"
puts $meshFile "Coordinates"
puts $meshFile "#node_number coord_x coord_y"
puts $meshFile " 1 0.000 0.000"
puts $meshFile " 2 2.000 0.000"
puts $meshFile " 3 2.000 2.000"
puts $meshFile " 4 0.000 2.000"
puts $meshFile " 10 2 4"
puts $meshFile " 11 0 4"
puts $meshFile " 16 4 0"
puts $meshFile " 17 4 2"
puts $meshFile " 22 4 4"
puts $meshFile " 26 6 0"
puts $meshFile " 27 6 2"
puts $meshFile " 32 6 4"
puts $meshFile " 36 8 0"
puts $meshFile " 37 8 2"
puts $meshFile " 42 8 4"

puts $meshFile "end coordinates"

puts $meshFile "Elements"
puts $meshFile "# element node1 node2 node3 node4"
puts $meshFile " 1 1 2 3 4"

puts $meshFile "2 4 3 10 11 "

puts $meshFile "3 2 16 17 3 "

puts $meshFile "4 3 17 22 10 "

puts $meshFile "5 16 26 27 17 "

puts $meshFile "6 17 27 32 22 "

puts $meshFile "7 26 36 37 27 "

puts $meshFile "8 27 37 42 32 "

puts $meshFile "end elements"

close $meshFile

set eleFile [open elementInfo.dat w]

puts $eleFile " 1 1 2 3 4"

close $eleFile

#-----------------------------------------------------------------------------------------
# 9. DEFINE ANALYSIS PARAMETERS
#-----------------------------------------------------------------------------------------

#---GROUND MOTION PARAMETERS
# time step in ground motion record

set motionDT 0.01

# number of steps in ground motion record

set motionSteps 2987

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

# element size (m)

set eleSize 2
# duration of ground motion (s)
set duration [expr $motionDT*$motionSteps]
# trial analysis time step
set kTrial [expr $eleSize/(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"

#-----------------------------------------------------------------------------------------
# 10. GRAVITY ANALYSIS
#-----------------------------------------------------------------------------------------

model BasicBuilder -ndm 2 -ndf 3

# loading object
pattern Plain 20 Constant {

load 11 0 -10 0

load 10 0 -20 0

load 22 0 -20 0

load 32 0 -20 0

load 42 0 -10 0

}

model BasicBuilder -ndm 2 -ndf 2

pattern Plain 30 Constant {

load 13 0 -40

load 24 0 -40

load 34 0 -40

load 44 0 -40

}

model BasicBuilder -ndm 2 -ndf 3

# update materials to ensure elastic behavior

updateMaterialStage -material 1 -stage 0

updateMaterialStage -material 2 -stage 0


set gamma 1.6

# create the SOE, ConstraintHandler, Integrator, Algorithm and Numberer

integrator Newmark $gamma [expr pow($gamma+0.5, 2)/4] 0.00 0.0 0.00 0.0

test EnergyIncr 1.0e-8 400 1;

constraints Transformation

algorithm Newton

numberer RCM

system ProfileSPD

analysis Transient

set startT [clock seconds]
analyze 10 5.0e2
analyze 10 5.0e3
puts "Finished with elastic gravity analysis..."

# update materials to consider plastic behavior

updateMaterialStage -material 1 -stage 1

updateMaterialStage -material 2 -stage 1

# plastic gravity loading
analyze 10 5.0e-3
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<=8} {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
updateParameter 10001 $vPerm2
updateParameter 10002 $hPerm2
updateParameter 10003 $vPerm2
updateParameter 10004 $hPerm2
updateParameter 10005 $vPerm1
updateParameter 10006 $hPerm1
updateParameter 10007 $vPerm1
updateParameter 10008 $hPerm1
updateParameter 10009 $vPerm1
updateParameter 10010 $hPerm1
updateParameter 10011 $vPerm1
updateParameter 10012 $hPerm1
updateParameter 10013 $vPerm2
updateParameter 10014 $hPerm2
updateParameter 10015 $vPerm2
updateParameter 10016 $hPerm2

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

recorder Element -file stress1.out -time -dT $recDT -eleRange 1 2 material 1 stress
recorder Element -file stress2.out -time -dT $recDT -eleRange 1 2 material 2 stress
recorder Element -file stress3.out -time -dT $recDT -eleRange 1 2 material 3 stress
recorder Element -file stress4.out -time -dT $recDT -eleRange 1 2 material 4 stress
recorder Element -file stress5.out -time -dT $recDT -eleRange 1 2 material 5 stress
recorder Element -file stress6.out -time -dT $recDT -eleRange 1 2 material 6 stress
recorder Element -file stress7.out -time -dT $recDT -eleRange 1 2 material 7 stress
recorder Element -file stress8.out -time -dT $recDT -eleRange 1 2 material 8 stress
recorder Element -file stress9.out -time -dT $recDT -eleRange 1 2 material 9 stress
recorder Element -file strain1.out -time -dT $recDT -eleRange 1 2 material 1 strain
recorder Element -file strain2.out -time -dT $recDT -eleRange 1 2 material 2 strain
recorder Element -file strain3.out -time -dT $recDT -eleRange 1 2 material 3 strain
recorder Element -file strain4.out -time -dT $recDT -eleRange 1 2 material 4 strain
recorder Element -file strain5.out -time -dT $recDT -eleRange 1 2 material 5 strain
recorder Element -file strain6.out -time -dT $recDT -eleRange 1 2 material 6 strain
recorder Element -file strain7.out -time -dT $recDT -eleRange 1 2 material 7 strain
recorder Element -file strain8.out -time -dT $recDT -eleRange 1 2 material 8 strain
recorder Element -file strain9.out -time -dT $recDT -eleRange 1 2 material 9 strain
puts "Finished creating all recorders..."

#-----------------------------------------------------------------------------------------
# 13. DYNAMIC ANALYSIS
#-----------------------------------------------------------------------------------------

model BasicBuilder -ndm 2 -ndf 3

# define constant scaling factor for applied velocity
set cFactor [expr $baseArea*$dashpotCoeff]

# define velocity time history file
set velocityFile velocityHistory.out

# timeseries object for force history
set mSeries "Path -dt $motionDT -filePath $velocityFile -factor $cFactor"

# loading object
pattern Plain 10 $mSeries {
load 1 1.0 0.0 0.0
}
puts "Dynamic loading created..."

set gamma 1.6

# create the SOE, ConstraintHandler, Integrator, Algorithm and Numberer

integrator Newmark $gamma [expr pow($gamma+0.5, 2)/4] 0.00 0.0 0.00 0.0

#test EnergyIncr 1.0e-8 400 1;

test NormDispIncr 1.0e-3 400 1

constraints Transformation

algorithm Newton

numberer RCM

system ProfileSPD

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]
puts "curTime: $curTime"
set curStep [expr $curTime/$dT]
puts "curStep: $curStep"
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/$dT]
puts "curStep: $curStep"
set remStep [expr int(($remStep-$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
Top
binger951
Posts: 1
Joined: Wed Aug 10, 2011 10:47 pm
Contact:
Contact binger951

复件1 15

Post by binger951 »

Lucky you. The first time I complained to them they dicked me around for six months and said they don't handle the type of problem I had (an unauthorized charge on my card).
lv sunglasses, chanel sunglasses, gucci sunglasses
Top
Locked

Return to “Soil Modelling”