Dynamic analysis with a modeled reservoir

[yingying]

I have meshed a dam and a reservoir (see image: https://ibb.co/nq0207). I am trying to do a coupled seismic analysis and I currently have the gravity analysis working with the following boundary conditions:

Base of dam is fixed in the vertical direction. Sides of foundation (orange in the image) are tied DOFs. The side of the reservoir (yellow) is fixed in the horizontal direction.

The results of the gravity analysis (displacement and pwp) look reasonable, however I cannot get a dynamic analysis to work.

I have tried changing the boundary conditions of the reservoir to Lysmer dashpots (this causes the gravity analysis to not run at all). I have a similar model, but without the reservoir, and a dynamic analysis works fine on it. The reason I have included the reservoir is so that the pwp and the hydrostatic pressures can be defined.

I have defined the reservoir using the following material "nDMaterial ElasticIsotropic 5 30 0.499998 1" and have created it using quadUP elements. Is this causing my dynamic analysis to not run at all?

My script is as follows:

#-------------------------------------------------------------------------------------------
# DEFINE NODES FOR SOIL ELEMENTS
#-------------------------------------------------------------------------------------------

# soil nodes are created in 2 dimensions, with 2 translational dof
model BasicBuilder -ndm 2 -ndf 3

# 2 ndm (number of dimensions), therefore NODE ID X Y
# define soil nodes

source nodes.tcl

puts "Finished creating all soil nodes..."


#-------------------------------------------------------------------------------------------
# DEFINE BOUNDARY CONDITIONS
#-------------------------------------------------------------------------------------------

# define fixity of base nodes for gravity analysis

source fix.tcl

# creating equal DOF for foundation layer base and sides
equalDOF 286 318 1 2
equalDOF 319 351 1 2
equalDOF 352 384 1 2
equalDOF 385 417 1 2
equalDOF 418 450 1 2
equalDOF 451 483 1 2
equalDOF 484 516 1 2
equalDOF 517 549 1 2
equalDOF 550 582 1 2
equalDOF 583 615 1 2
equalDOF 616 648 1 2
equalDOF 649 681 1 2
equalDOF 682 714 1 2
equalDOF 715 747 1 2
equalDOF 748 780 1 2
equalDOF 781 813 1 2
equalDOF 814 846 1 2
equalDOF 847 879 1 2
equalDOF 880 912 1 2
equalDOF 913 945 1 2
equalDOF 946 963 1 2


# pwp

fix 218 0 0 1
fix 219 0 0 1
fix 220 0 0 1
fix 221 0 0 1
fix 83 0 0 1
fix 58 0 0 1
fix 59 0 0 1
fix 60 0 0 1
fix 104 0 0 1
fix 278 0 0 1
fix 279 0 0 1
fix 280 0 0 1
fix 281 0 0 1
fix 222 0 0 1
fix 223 0 0 1
fix 224 0 0 1
fix 84 0 0 1
fix 61 0 0 1
fix 63 0 0 1
fix 105 0 0 1
fix 282 0 0 1
fix 283 0 0 1
fix 284 0 0 1
fix 285 0 0 1

fix 150 0 0 1
fix 155 0 0 1
fix 160 0 0 1
fix 165 0 0 1
fix 164 0 0 1
fix 229 0 0 1
fix 233 0 0 1
fix 237 0 0 1
fix 241 0 0 1
fix 245 0 0 1
fix 249 0 0 1
fix 253 0 0 1
fix 257 0 0 1
fix 261 0 0 1
fix 265 0 0 1
fix 269 0 0 1
fix 273 0 0 1
fix 277 0 0 1

fix 145 0 0 1
fix 140 0 0 1
fix 955 0 0 1
fix 956 0 0 1
fix 957 0 0 1
fix 958 0 0 1
fix 959 0 0 1
fix 960 0 0 1
fix 961 0 0 1
fix 962 0 0 1
fix 963 0 0 1


puts "Finished creating all boundary conditions..."

#-------------------------------------------------------------------------------------------
# MESH GEOMETRY
#-------------------------------------------------------------------------------------------

# number of nodes
set numNodes 1138

# number of elements
set numEle 1065

#-------------------------------------------------------------------------------------------
# MATERIAL PARAMETERS
#-------------------------------------------------------------------------------------------

# combined undrained bulk modulus (kPa)
set bulk [expr 2.2e6/1]

# fluid mass density
set fmass 1

# initial permeability coefficient in horizontal direction for gravity analysis
set hPerm 1

# initial permeability coefficient in vertical direction for gravity analysis
set vPerm 1

#peak shear strain (gamma max)
set peakSS 0.1

#Reference mean effective confining pressure (default is 100) (p'r)
set refPress 1200

#pressure dependent coefficient (default 0) (d)
set refCoef 2

# phase transformation angle (deg)

set phaseAng 27

#number of yield surfaces (default 20)
set noYieldSurf 25

#-------------------------------------------------------------------------------------------
# OTHER PROPERTIES
#-------------------------------------------------------------------------------------------

# USER DEFINED VARIABLES

# Gravity m/s2
set g 9.81

set gravX 0

set gravY [expr -$g]

#---RAYLEIGH DAMPING PARAMETERS
set pi 3.141592654
# damping ratio - 5 percent damping
set damp 0.05
# 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

#Analysis settings
#time step
set dT 0.02
#number of steps
set nSteps 400

#-------------------------------------------------------------------------------------------
# DEFINE MATERIAL
#-------------------------------------------------------------------------------------------


# clay core

set hPerm1 1e-10

set vPerm1 1e-10

#saturated soil mass density
set rho1 2

# shear modulus - 200x3.5 (kPa)
set G1 700000

#bulk modulus - roughly triple shear modulus (kPa)
set B1 1517000

# friction angle
set frictAng1 37

nDMaterial PressureDependMultiYield 1 2 $rho1 $G1 $B1 $frictAng1 $peakSS $refPress $refCoef $phaseAng 0.0 0.0 0.0 0.0 0.0 0.0 $noYieldSurf




# sand filters

set hPerm2 1e-3

set vPerm2 1e-3

#saturated soil mass density
set rho2 2.18

# shear modulus - 200x3.5 (kPa)
set G2 700000

#bulk modulus - roughly triple shear modulus (kPa)
set B2 1517000

# friction angle
set frictAng2 37

nDMaterial PressureDependMultiYield 2 2 $rho2 $G2 $B2 $frictAng2 $peakSS $refPress $refCoef $phaseAng 0.0 0.0 0.0 0.0 0.0 0.0 $noYieldSurf





#rockfill shells

set hPerm3 1e-2

set vPerm3 1e-2

#saturated soil mass density
set rho3 2.08

# shear modulus - 200x3.5 (kPa)
set G3 700000

#bulk modulus - roughly triple shear modulus (kPa)
set B3 1517000

# friction angle
set frictAng3 37

nDMaterial PressureDependMultiYield 3 2 $rho3 $G3 $B3 $frictAng3 $peakSS $refPress $refCoef $phaseAng 0.0 0.0 0.0 0.0 0.0 0.0 $noYieldSurf




# foundation alluvium

set hPerm4 1e-7

set vPerm4 1e-7

#saturated soil mass density
set rho4 2.08

# shear modulus - 200x3.5 (kPa)
set G4 700000

#bulk modulus - roughly triple shear modulus (kPa)
set B4 1517000

# friction angle
set frictAng4 37

nDMaterial PressureDependMultiYield 4 2 $rho4 $G4 $B4 $frictAng4 $peakSS $refPress $refCoef $phaseAng 0.0 0.0 0.0 0.0 0.0 0.0 $noYieldSurf


# reservoir water

nDMaterial ElasticIsotropic 5 30 0.499998 1

puts "Finished creating materials..."

#-------------------------------------------------------------------------------------------
# DASHPOT
#-------------------------------------------------------------------------------------------

source dashpot.tcl

# dashpot size element x rock density x rock VS
set mC 1000

uniaxialMaterial Viscous 4000 $mC 1

element zeroLength 3000 2000 2001 -mat 4000 -dir 1
element zeroLength 3001 2002 2003 -mat 4000 -dir 1
element zeroLength 3002 2004 2005 -mat 4000 -dir 1
element zeroLength 3003 2006 2007 -mat 4000 -dir 1
element zeroLength 3004 2008 2009 -mat 4000 -dir 1
element zeroLength 3005 2010 2011 -mat 4000 -dir 1
element zeroLength 3006 2012 2013 -mat 4000 -dir 1
element zeroLength 3007 2014 2015 -mat 4000 -dir 1
element zeroLength 3008 2016 2017 -mat 4000 -dir 1
element zeroLength 3009 2018 2019 -mat 4000 -dir 1
element zeroLength 3010 2020 2021 -mat 4000 -dir 1
element zeroLength 3011 2022 2023 -mat 4000 -dir 1
element zeroLength 3012 2024 2025 -mat 4000 -dir 1
element zeroLength 3013 2026 2027 -mat 4000 -dir 1
element zeroLength 3014 2028 2029 -mat 4000 -dir 1
element zeroLength 3015 2030 2031 -mat 4000 -dir 1
element zeroLength 3016 2032 2033 -mat 4000 -dir 1
element zeroLength 3017 2034 2035 -mat 4000 -dir 1

puts "Finished creating dashpot material..."



#-------------------------------------------------------------------------------------------
# DEFINE SOIL ELEMENTS
#-------------------------------------------------------------------------------------------

# define soil elements

source ele.tcl

puts "Finished creating all elements..."


#------------------------------------------------------------------------------------------
# CREATE GRAVITY RECORDERS
#------------------------------------------------------------------------------------------

# record nodal displacement and acceleration

recorder Node -file Gdisp.out -time -nodeRange 1 $numNodes -dof 1 2 disp
recorder Node -file Gpwp.out -time -nodeRange 1 $numNodes -dof 3 vel

puts "Finished creating gravity recorders..."


#-------------------------------------------------------------------------------------------
# APPLY GRAVITY LOADING
#-------------------------------------------------------------------------------------------

# update materials to ensure elastic behavior
updateMaterialStage -material 1 -stage 0
updateMaterialStage -material 2 -stage 0
updateMaterialStage -material 3 -stage 0
updateMaterialStage -material 4 -stage 0


timeSeries Linear 1499

pattern Plain 3 1499 {
sp 946 3 529.74
sp 947 3 529.74
sp 948 3 529.74
sp 949 3 529.74
sp 950 3 529.74
sp 951 3 529.74
sp 952 3 529.74
sp 953 3 529.74
sp 954 3 529.74
sp 106 3 529.74
sp 111 3 500.31
sp 116 3 470.88
sp 121 3 441.45
sp 126 3 412.02
sp 131 3 382.59
sp 132 3 382.59
sp 166 3 353.16
sp 170 3 323.73
sp 174 3 294.3
sp 178 3 264.87
sp 182 3 235.44
sp 186 3 206.01
sp 190 3 176.58
sp 194 3 147.15
sp 198 3 117.72
sp 202 3 88.29
sp 206 3 58.86
sp 210 3 29.43
sp 214 3 0

}


constraints Transformation
test NormDispIncr 0.0005 30 2
algorithm Newton
numberer RCM
system ProfileSPD
integrator Newmark 1.5 1
analysis Transient

analyze 10 100

puts "Finished with elastic gravity analysis..."


# update material to consider plastic behavior
updateMaterialStage -material 1 -stage 1
updateMaterialStage -material 2 -stage 1
updateMaterialStage -material 3 -stage 1
updateMaterialStage -material 4 -stage 1


# execute plastic gravity analysis
analyze 10 0.005

# holds gravity constant and restart time

puts "Finished with plastic gravity analysis..."


#------------------------------------------------------------------------------------------
# UPDATE ELEMENT PERMEABILITY
#-------------------------------------------------------------------------------------------

source Paramh_r_coup.tcl

source Paramv_r_coup.tcl


source updateParamh_r_coup.tcl

puts "Finished updating permh..."


source updateParamv_r_coup.tcl

puts "Finished updating permv..."


#------------------------------------------------------------------------------------------
# CREATE DYNAMIC RECORDERS
#-------------------------------------------------------------------------------------------

# reset time and analysis
setTime 0.0
wipeAnalysis
remove recorders

# want to record acceleration
recorder Node -file disp.out -time -nodeRange 1 $numNodes -dof 1 2 disp
#recorder Node -file vel.out -time -nodeRange 1 $numNodes -dof 1 2 vel
#recorder Node -file accelbase.out -time -node 1071 -dof 1 2 accel
recorder Node -file accelcrest.out -time -node 62 -dof 1 accel
recorder Node -file pwp.out -time -nodeRange 1 $numNodes -dof 3 vel


#recorder Element -file stress1.out -time -ele 103 material 1 stress
#recorder Element -file strain1.out -time -ele 103 material 1 strain


#-------------------------------------------------------------------------------------------
# APPLY EARTHQUAKE ACCELERATION
#------------------------------------------------------------------------------------------

#acceleration file
set accelFile "eq2accel.csv"
set cFactor 0.00001

#timeseries object
#set accelSeries "Series -dt $dT -filePath $accelFile -factor $cFactor"

timeSeries Path 1500 -dt $dT -filePath $accelFile -factor $cFactor

#apply acceleration load (tagID) (direction of excitation)
pattern UniformExcitation 2 1 -accel 1500

puts "Dynamic loading created..."

# clear previously defined analysis parameters
wipeAnalysis

constraints Penalty 1.e18 1.e18
test NormDispIncr 0.005 50 2
algorithm KrylovNewton
numberer RCM
system ProfileSPD
rayleigh $a0 $a1 0.0 0.0
integrator Newmark $gamma $beta
analysis Transient

# executes the analysis
analyze $nSteps $dT

puts "Finished with plastic dynamic analysis..."


wipe

[jawaheer]

HI
were you able to solve the problem?
Thanks for your answer