A problem about mass density of quad and brick element

Forum for OpenSees users to post questions, comments, etc. on the use of the OpenSees interpreter, OpenSees.exe

Moderators: silvia, selimgunay, Moderators

Post Reply
nvzitejing
Posts: 6
Joined: Wed Aug 15, 2012 11:27 pm
Location: Tongji University

A problem about mass density of quad and brick element

Post by nvzitejing »

Hi,everyone!
I am analyzing site response problem in OpenSees 2.4.0. Referring to the examples"2D Total Stress Site Response Analysis of a Layered Soil Column"and "3D Site Response Analysis of Sloping Ground", I mode layers soil using nDmaterial and quad or brick element. In the nDmaterial and some quad or brick element models, the mass density will be included. But when I perform eigen analysis for the model, the period results of 3D are 1.414 times ones of 2D and the ones obtained using ANSYS software. I set mass density as 0.0 and use node mass for model, then perform eigen analysis again. The results of 3D and 2D are the same as ones of ANSYS.
I don't understand why is that so ? What's the matter with distributed mass? would you please help me?


2D model as following:

wipe

set dataDir Data-2ndm; # set up name of data directory -- remove
set dataDir1 SoilData;
set dataDirSoil $dataDir/$dataDir1; # set up name of data directory -- remove
file mkdir $dataDirSoil;
source LibUnits.tcl;


#-------------------------------------------------------------------------------------------
#1.SOIL GEOMETRY
#-------------------------------------------------------------------------------------------
#thickness of soil profile (m)
set soilThick 116.8
# number of soil layers
set numLayers 21
# layer thicknessess (m)
set layerThick(21) 7.0
set layerThick(20) 5.0
set layerThick(19) 7.2
set layerThick(18) 7.4
set layerThick(17) 5.0
set layerThick(16) 5.0
set layerThick(15) 9.9
set layerThick(14) 2.7
set layerThick(13) 6.0
set layerThick(12) 6.8
set layerThick(11) 2.7
set layerThick(10) 6.0
set layerThick(9) 4.6
set layerThick(8) 2.7
set layerThick(7) 8.4
set layerThick(6) 3.0
set layerThick(5) 6.6
set layerThick(4) 6.6
set layerThick(3) 4.0
set layerThick(2) 5.0
set layerThick(1) 5.2
#-------------------------------------------------------------------------------------------
#2.MESH GEOMETRY
#-------------------------------------------------------------------------------------------
# number of elements in horizontal direction
set nElemX 1
# horizontal element size (m)
set sElemX 10
# number of elements in vertical direction for each layer
set nElemY(21) 7
set nElemY(20) 5
set nElemY(19) 8
set nElemY(18) 8
set nElemY(17) 5
set nElemY(16) 5
set nElemY(15) 10
set nElemY(14) 3
set nElemY(13) 6
set nElemY(12) 8
set nElemY(11) 3
set nElemY(10) 5
set nElemY(9) 4
set nElemY(8) 3
set nElemY(7) 8
set nElemY(6) 3
set nElemY(5) 6
set nElemY(4) 6
set nElemY(3) 4
set nElemY(2) 5
set nElemY(1) 5
# total number 10f elements in vertical direction
set nElemT 0
for {set i 1} {$i<=$numLayers} {incr i 1} {
set nElemT [expr $nElemT+$nElemY($i)]
}
puts "$nElemT"
# vertical element size in each layer
for {set i 1} {$i <=$numLayers} {incr i 1} {
set sElemY($i) [expr $layerThick($i)/$nElemY($i)]
puts "size: $sElemY($i)"
}
# number of nodes in vertical direction in each layer
set nNodeT 0
for {set k 1} {$k<$numLayers} {incr k 1} {
set nNodeY($k) [expr 2*$nElemY($k)]
puts "number of nodes in layer $k: $nNodeY($k)"
set nNodeT [expr $nNodeT+$nNodeY($k)]
}
set nNodeY($numLayers) [expr 2*($nElemY($numLayers)+1)]
puts "number of nodes in layer $numLayers: $nNodeY($numLayers)"
set nNodeT [expr $nNodeT + $nNodeY($numLayers)]
puts "total number of nodes: $nNodeT"
#-----------------------------------------------------------------------------------------
#3.CREATE SOIL NODES
#-----------------------------------------------------------------------------------------
model BasicBuilder -ndm 2 -ndf 2
set yCoord 0.0
set count 0
set fileId1 [open $dataDirSoil/SoilNode.txt w ]
set fileId2 [open $dataDirSoil/yCoord.txt w ]
for {set k 1} {$k<=$numLayers} {incr k 1} {
for {set j 1} {$j<=$nNodeY($k)} {incr j 2} {
node [expr $j+$count] 0.0 $yCoord
node [expr $j+$count+1] $sElemX $yCoord


puts $fileId1 " node [expr $j+$count] 0.0 $yCoord "
puts $fileId1 " node [expr $j+$count+1] $sElemX $yCoord "
set yCoord [expr $yCoord+$sElemY($k)]
puts $fileId2 " node [expr $j+$count] $yCoord "
}
set count [expr $count+$nNodeY($k)]
}
close $fileId1
close $fileId2
puts "Finished creating all soil nodes"

# define boundary conditions for nodes at base of column
fix 1 1 1
fix 2 1 1

# define periodic boundary conditions for remaining nodes
for {set k 3} {$k <=[expr $nNodeT]} {incr k 2} {
equalDOF $k [expr $k+1] 1 2
puts "equalDOF $k [expr $k+1] 1 2 "
}
puts "Finished creating all soil boundary conditions..."
#-----------------------------------------------------------------------------------------
#4.CREATE SOIL MATERIALS
#-----------------------------------------------------------------------------------------
#---MATERIAL PROPERTIES
# soil mass density (Mg/m^3)
set rho(21) 1.95
set rho(20) 2.01
set rho(19) 2.00
set rho(18) 2.01
set rho(17) 2.00
set rho(16) 2.09
set rho(15) 2.16
set rho(14) 2.16
set rho(13) 2.11
set rho(12) 2.11
set rho(11) 2.11
set rho(10) 2.11
set rho(9) 2.28
set rho(8) 2.28
set rho(7) 2.28
set rho(6) 2.28
set rho(5) 2.02
set rho(4) 2.16
set rho(3) 2.07
set rho(2) 2.07
set rho(1) 2.07
# soil shear wave velocity (m/s)
set Vs(21) 136
set Vs(20) 145
set Vs(19) 147
set Vs(18) 160
set Vs(17) 175
set Vs(16) 201
set Vs(15) 245
set Vs(14) 315
set Vs(13) 301
set Vs(12) 341
set Vs(11) 357
set Vs(10) 349
set Vs(9) 324
set Vs(8) 347
set Vs(7) 338
set Vs(6) 341
set Vs(5) 332
set Vs(4) 352
set Vs(3) 414
set Vs(2) 459
set Vs(1) 534
# soil shear modulus (kPa)
for {set i 1} {$i<=$numLayers} {incr i 1} {
set G($i) [expr $rho($i)*$Vs($i)*$Vs($i)]
}
# poisson's ratio of soil
set nu(21) 0.471
set nu(20) 0.469
set nu(19) 0.471
set nu(18) 0.461
set nu(17) 0.448
set nu(16) 0.436
set nu(15) 0.423
set nu(14) 0.381
set nu(13) 0.394
set nu(12) 0.362
set nu(11) 0.39
set nu(10) 0.359
set nu(9) 0.382
set nu(8) 0.364
set nu(7) 0.378
set nu(6) 0.363
set nu(5) 0.372
set nu(4) 0.377
set nu(3) 0.351
set nu(2) 0.318
set nu(1) 0.33
# soil elastic modulus (kPa)
for {set i 1} {$i<=$numLayers} {incr i 1} {
set E($i) [expr 2*$G($i)*(1+$nu($i))]
}
# soil bulk modulus (kPa)
for {set i 1} {$i<=$numLayers} {incr i 1} {
set Bulk($i) [expr $E($i)/(3*(1-2*$nu($i)))]
}
# soil cohesion (kPa)
set Cohesion(5) 73
set Cohesion(2) 90
set Cohesion(1) 115
# peak shear strain
set gammaPeak 0.1
# soil friction angle
set phi(21) 32.0
set phi(20) 34.2
set phi(19) 34.2
set phi(18) 28.4
set phi(17) 34.1
set phi(16) 34.7
set phi(15) 37.8
set phi(14) 30.3
set phi(13) 36.1
set phi(12) 38.5
set phi(11) 35.0
set phi(10) 36.6
set phi(9) 36.6
set phi(8) 36.7
set phi(7) 36.0
set phi(6) 35.0
set phi(5) 0.0
set phi(4) 33.0
set phi(3) 34.3
set phi(2) 0.0
set phi(1) 0.0
#phase transformation angle
set phaseAng 27
# reference pressure
set refPress 80.0
# pressure dependency coefficient
set pressCoeff 0
# contraction
set contract 0.06
# dilation coefficients
set dilate1 0.6
set dilate2 3
# liquefaction coefficients
set liq1 0
set liq2 0
set liq3 0
# bedrock shear wave velocity (m/s)
set rockVS 600
# bedrock mass density (Mg/m^3)
set rockDen 2.4
nDMaterial PressureIndependMultiYield 1 2 0.0 $G(1) $Bulk(1) $Cohesion(1) $gammaPeak \
$phi(1) $refPress $pressCoeff 30
nDMaterial PressureIndependMultiYield 2 2 0.0 $G(2) $Bulk(2) $Cohesion(2) $gammaPeak \
$phi(2) $refPress $pressCoeff 30

for {set k 3} {$k<=4} {incr k 1} {
nDMaterial PressureDependMultiYield $k 2 0.0 $G($k) $Bulk($k) $phi($k) $gammaPeak $refPress $pressCoeff\
$phaseAng $contract $dilate1 $dilate2 $liq1 $liq2 $liq3 30
}
nDMaterial PressureIndependMultiYield 5 2 0.0 $G(5) $Bulk(5) $Cohesion(5) $gammaPeak \
$phi(5) $refPress $pressCoeff 30
for {set k 6} {$k<=21} {incr k 1} {
nDMaterial PressureDependMultiYield $k 2 0.0 $G($k) $Bulk($k) $phi($k) $gammaPeak $refPress $pressCoeff\
$phaseAng $contract $dilate1 $dilate2 $liq1 $liq2 $liq3 30
}
#for {set k 1} {$k<=21} {incr k 1} {
# nDMaterial ElasticIsotropic $k $E($k) $nu($k) 0.0
#}
puts "Finished soil material define......"
#-----------------------------------------------------------------------------------------
#5.NODE MASS OF SOIL
#-----------------------------------------------------------------------------------------
#Elemass----------------------------------------------------------------
set count 0
set fileId [open $dataDirSoil/Elemass.txt w ]
for {set k 1} {$k<=$numLayers} {incr k 1} {
for {set j 1} {$j<=$nElemY($k)} {incr j 1} {
set Elemass([expr $j+$count]) [expr $rho($k)*$sElemY($k)*$sElemX]
puts $fileId " [expr $j+$count] $Elemass([expr $j+$count]) "
}
set count [expr $count+$nElemY($k)]
}
close $fileId
set fileId [open $dataDirSoil/Nodemass.txt w ]
mass 1 [expr $Elemass(1)/4] [expr $Elemass(1)/4]
mass 2 [expr $Elemass(1)/4] [expr $Elemass(1)/4]
puts $fileId "node 1-2 [expr $Elemass(1)/4] "
for {set i 1} {$i<=[expr $nElemT-1] } {incr i 1} {
mass [expr 2*$i+1] [expr ($Elemass($i)+$Elemass([expr $i+1]))/4] [expr ($Elemass($i)+$Elemass([expr $i+1]))/4]
mass [expr 2*$i+2] [expr ($Elemass($i)+$Elemass([expr $i+1]))/4] [expr ($Elemass($i)+$Elemass([expr $i+1]))/4]
puts $fileId " node [expr 2*$i+1]-[expr 2*$i+2] [expr ($Elemass($i)+$Elemass([expr $i+1]))/4] "
}
mass [expr 2*$nElemT+1] [expr $Elemass($nElemT)/4] [expr $Elemass($nElemT)/4]
mass [expr 2*$nElemT+2] [expr $Elemass($nElemT)/4] [expr $Elemass($nElemT)/4]
puts $fileId "node [expr 2*$nElemT+1]-[expr 2*$nElemT+2] [expr $Elemass(1)/4]"
close $fileId
#-----------------------------------------------------------------------------------------
#6.CREATE SOIL ELEMENTS
#-----------------------------------------------------------------------------------------

set wgtX 0.0
for {set k 1} {$k<=$numLayers} {incr k 1} {
set wgtY($k) [expr -9.81*$rho($k)]
}
set count 0
set fileId [open $dataDirSoil/SoilElements.txt w ]
# loop over layers
for {set k 1} {$k<=$numLayers} {incr k 1} {
# loop over elements
for {set j 1} {$j<=$nElemY($k)} {incr j 1} {
set nI [expr 2*($j+$count)-1]
set nJ [expr $nI+1]
set nK [expr $nI+3]
set nL [expr $nI+2]
element quad [expr $j+$count] $nI $nJ $nK $nL 1.0 "PlaneStrain" $k 0.0 0.0 $wgtX $wgtY($k)
puts $fileId "[expr $j+$count] $nI $nJ $nK $nL"
}
set count [expr $count+$nElemY($k)]
}
close $fileId
puts "Finished creating all soil elements..."

set numModes 10;
for { set k 1 } { $k<=$numModes } { incr k 1 } {
recorder Node -file $dataDirSoil/mode$k-1.out -nodeRange 1 236 -dof 1 "eigen $k"
recorder Node -file $dataDirSoil/mode$k-2.out -nodeRange 1 236 -dof 2 "eigen $k"
}

set lambda [eigen $numModes];

set omega {};
set Freq {};
set T {};

foreach lam $lambda {
lappend omega [expr sqrt($lam)]
lappend Freq [expr sqrt($lam)/(2*$pi)]
lappend T [expr (2*$pi)/sqrt($lam)]
}

# Outputting Freq. & Period
set Frequency [open $dataDirSoil/Freq.txt "w"]
foreach t $Freq {
puts $Frequency " $t"
}
close $Frequency

set Periods [open $dataDirSoil/period.txt "w"]
foreach t $T {
puts $Periods " $t"
}
close $Periods
puts " Analysis of EigenValue Done"



And 3D model as following :
wipe

set dataDir Data; # set up name of data directory -- remove
set dataDir1 SoilData;
set dataDirSoil $dataDir/$dataDir1; # set up name of data directory -- remove
file mkdir $dataDirSoil;
source LibUnits.tcl;

#-------------------------------------------------------------------------------------------
#1.SOIL GEOMETRY
#-------------------------------------------------------------------------------------------
#thickness of soil profile (m)
set soilThick 116.8
# number of soil layers
set numLayers 21
# layer thicknessess (m)
set layerThick(21) 7.0
set layerThick(20) 5.0
set layerThick(19) 7.2
set layerThick(18) 7.4
set layerThick(17) 5.0
set layerThick(16) 5.0
set layerThick(15) 9.9
set layerThick(14) 2.7
set layerThick(13) 6.0
set layerThick(12) 6.8
set layerThick(11) 2.7
set layerThick(10) 6.0
set layerThick(9) 4.6
set layerThick(8) 2.7
set layerThick(7) 8.4
set layerThick(6) 3.0
set layerThick(5) 6.6
set layerThick(4) 6.6
set layerThick(3) 4.0
set layerThick(2) 5.0
set layerThick(1) 5.2
# define layer boundaries
set layerBound(1) $layerThick(1)
for {set i 2} {$i<=$numLayers} {incr i 1} {
set layerBound($i) [expr $layerBound([expr $i-1])+$layerThick($i)]
}
#-------------------------------------------------------------------------------------------
#2.MESH GEOMETRY
#-------------------------------------------------------------------------------------------
# number of elements in horizontal direction
set nElemX 1
set nElemZ 1
# number of nodes in horizontal direction
set nNodeX [expr 2*$nElemX+1]
# horizontal element size (m)
set sElemX 1.0
set sElemZ 1.0
# number of elements in vertical direction for each layer
set nElemY(21) 7
set nElemY(20) 5
set nElemY(19) 8
set nElemY(18) 8
set nElemY(17) 5
set nElemY(16) 5
set nElemY(15) 10
set nElemY(14) 3
set nElemY(13) 6
set nElemY(12) 8
set nElemY(11) 3
set nElemY(10) 5
set nElemY(9) 4
set nElemY(8) 3
set nElemY(7) 8
set nElemY(6) 3
set nElemY(5) 6
set nElemY(4) 6
set nElemY(3) 4
set nElemY(2) 5
set nElemY(1) 5
# total number of elements in vertical direction
set nElemT 0
for {set i 1} {$i<=$numLayers} {incr i 1} {
set nElemT [expr $nElemT+$nElemY($i)]
}
puts "$nElemT"
# vertical element size in each layer
for {set i 1} {$i <=$numLayers} {incr i 1} {
set sElemY($i) [expr $layerThick($i)/$nElemY($i)]
puts "size: $sElemY($i)"
}
# number of nodes in vertical direction in each layer
set nNodeT 0
for {set k 1} {$k<$numLayers} {incr k 1} {
set nNodeY($k) [expr 4*$nElemY($k)]
puts "number of nodes in layer $k: $nNodeY($k)"
set nNodeT [expr $nNodeT+$nNodeY($k)]
}
set nNodeY($numLayers) [expr 4*($nElemY($numLayers)+1)]
puts "number of nodes in layer $numLayers: $nNodeY($numLayers)"
set nNodeT [expr $nNodeT + $nNodeY($numLayers)]
puts "total number of nodes: $nNodeT"
#-----------------------------------------------------------------------------------------
#3.CREATE SOIL NODES
#-----------------------------------------------------------------------------------------
model BasicBuilder -ndm 3 -ndf 3
set yCoord 0.0
set count 0
set fileId [open $dataDirSoil/SoilNode.txt w ]
for {set k 1} {$k<=$numLayers} {incr k 1} {
for {set j 1} {$j<=$nNodeY($k)} {incr j 4} {
node [expr $j+$count] $sElemX $yCoord $sElemZ
node [expr $j+$count+1] $sElemX $yCoord 0.0
node [expr $j+$count+2] 0.0 $yCoord 0.0
node [expr $j+$count+3] 0.0 $yCoord $sElemZ


puts $fileId " node [expr $j+$count] $sElemX $yCoord $sElemZ"
puts $fileId " node [expr $j+$count+1] $sElemX $yCoord 0.0"
puts $fileId " node [expr $j+$count+2] 0.0 $yCoord 0.0"
puts $fileId " node [expr $j+$count+3] 0.0 $yCoord $sElemZ"
set yCoord [expr $yCoord+$sElemY($k)]
}
set count [expr $count + $nNodeY($k)]
}
close $fileId
puts "Finished creating all soil nodes"

# define boundary conditions for nodes at base of column
fix 1 1 1 1
fix 2 1 1 1
fix 3 1 1 1
fix 4 1 1 1

# define periodic boundary conditions for remaining nodes
for {set k 5} {$k <= [expr $nNodeT]} {incr k 4} {
equalDOF $k [expr $k+1] 1 2 3
equalDOF $k [expr $k+2] 1 2 3
equalDOF $k [expr $k+3] 1 2 3
puts "equalDOF $k [expr $k+1] 1 2 3"
puts "equalDOF $k [expr $k+2] 1 2 3"
puts "equalDOF $k [expr $k+3] 1 2 3"
}
puts "Finished creating all soil boundary conditions..."
#-----------------------------------------------------------------------------------------
#4.CREATE SOIL MATERIALS
#-----------------------------------------------------------------------------------------
#---MATERIAL PROPERTIES
# soil mass density (Mg/m^3)
set rho(21) 1.95
set rho(20) 2.01
set rho(19) 2.00
set rho(18) 2.01
set rho(17) 2.00
set rho(16) 2.09
set rho(15) 2.16
set rho(14) 2.16
set rho(13) 2.11
set rho(12) 2.11
set rho(11) 2.11
set rho(10) 2.11
set rho(9) 2.28
set rho(8) 2.28
set rho(7) 2.28
set rho(6) 2.28
set rho(5) 2.02
set rho(4) 2.16
set rho(3) 2.07
set rho(2) 2.07
set rho(1) 2.07
# soil shear wave velocity (m/s)
set Vs(21) 136
set Vs(20) 145
set Vs(19) 147
set Vs(18) 160
set Vs(17) 175
set Vs(16) 201
set Vs(15) 245
set Vs(14) 315
set Vs(13) 301
set Vs(12) 341
set Vs(11) 357
set Vs(10) 349
set Vs(9) 324
set Vs(8) 347
set Vs(7) 338
set Vs(6) 341
set Vs(5) 332
set Vs(4) 352
set Vs(3) 414
set Vs(2) 459
set Vs(1) 534
# soil shear modulus (kPa)
for {set i 1} {$i<=$numLayers} {incr i 1} {
set G($i) [expr $rho($i)*$Vs($i)*$Vs($i)]
}
# poisson's ratio of soil
set nu(21) 0.471
set nu(20) 0.469
set nu(19) 0.471
set nu(18) 0.461
set nu(17) 0.448
set nu(16) 0.436
set nu(15) 0.423
set nu(14) 0.381
set nu(13) 0.394
set nu(12) 0.362
set nu(11) 0.39
set nu(10) 0.359
set nu(9) 0.382
set nu(8) 0.364
set nu(7) 0.378
set nu(6) 0.363
set nu(5) 0.372
set nu(4) 0.377
set nu(3) 0.351
set nu(2) 0.318
set nu(1) 0.33
#set mu 0.0
# soil elastic modulus (kPa)
for {set i 1} {$i<=$numLayers} {incr i 1} {
set E($i) [expr 2*$G($i)*(1+$nu($i))]
}
# soil bulk modulus (kPa)
for {set i 1} {$i<=$numLayers} {incr i 1} {
set Bulk($i) [expr $E($i)/(3*(1-2*$nu($i)))]
}
# soil cohesion (kPa)
set Cohesion(5) 73
set Cohesion(2) 90
set Cohesion(1) 115
# peak shear strain
set gammaPeak 0.1
# soil friction angle
set phi(21) 32.0
set phi(20) 34.2
set phi(19) 34.2
set phi(18) 28.4
set phi(17) 34.1
set phi(16) 34.7
set phi(15) 37.8
set phi(14) 30.3
set phi(13) 36.1
set phi(12) 38.5
set phi(11) 35.0
set phi(10) 36.6
set phi(9) 36.6
set phi(8) 36.7
set phi(7) 36.0
set phi(6) 35.0
set phi(5) 0.0
set phi(4) 33.0
set phi(3) 34.3
set phi(2) 0.0
set phi(1) 0.0
#phase transformation angle
set phaseAng 27
# reference pressure
set refPress 80.0
# pressure dependency coefficient
set pressCoeff 0
# contraction
set contract 0.06
# dilation coefficients
set dilate1 0.6
set dilate2 3
# liquefaction coefficients
set liq1 0
set liq2 0
set liq3 0
# bedrock shear wave velocity (m/s)
set rockVS 600
# bedrock mass density (Mg/m^3)
set rockDen 2.4

nDMaterial PressureIndependMultiYield 1 3 0 $G(1) $Bulk(1) $Cohesion(1) $gammaPeak \
$phi(1) $refPress $pressCoeff 30
nDMaterial PressureIndependMultiYield 2 3 0 $G(2) $Bulk(2) $Cohesion(2) $gammaPeak \
$phi(2) $refPress $pressCoeff 30

for {set k 3} {$k<=4} {incr k 1} {
nDMaterial PressureDependMultiYield $k 3 0 $G($k) $Bulk($k) $phi($k) $gammaPeak $refPress $pressCoeff\
$phaseAng $contract $dilate1 $dilate2 $liq1 $liq2 $liq3 30
}
nDMaterial PressureIndependMultiYield 5 3 0 $G(5) $Bulk(5) $Cohesion(5) $gammaPeak \
$phi(5) $refPress $pressCoeff 30
for {set k 6} {$k<=21} {incr k 1} {
nDMaterial PressureDependMultiYield $k 3 0 $G($k) $Bulk($k) $phi($k) $gammaPeak $refPress $pressCoeff\
$phaseAng $contract $dilate1 $dilate2 $liq1 $liq2 $liq3 30
}
#for {set k 1} {$k<=21} {incr k 1} {
# nDMaterial ElasticIsotropic $k $E($k) $nu($k) 0
#}
puts "Finished soil material define......"
#-----------------------------------------------------------------------------------------
#5.NODE MASS OF SOIL
#-----------------------------------------------------------------------------------------
#Elemass----------------------------------------------------------------
set count 0
set fileId [open $dataDirSoil/Elemass.txt w ]
for {set k 1} {$k<=$numLayers} {incr k 1} {
for {set j 1} {$j<=$nElemY($k)} {incr j 1} {
set Elemass([expr $j+$count]) [expr $rho($k)*$sElemY($k)*$sElemX]
puts $fileId " [expr $j+$count] $Elemass([expr $j+$count]) "
}
set count [expr $count+$nElemY($k)]
}
close $fileId
set fileId [open $dataDirSoil/Nodemass.txt w ]
mass 1 [expr $Elemass(1)/8] [expr $Elemass(1)/8] [expr $Elemass(1)/8]
mass 2 [expr $Elemass(1)/8] [expr $Elemass(1)/8] [expr $Elemass(1)/8]
mass 3 [expr $Elemass(1)/8] [expr $Elemass(1)/8] [expr $Elemass(1)/8]
mass 4 [expr $Elemass(1)/8] [expr $Elemass(1)/8] [expr $Elemass(1)/8]
puts $fileId "node 1-4 [expr $Elemass(1)/8] "
for {set i 1} {$i<=[expr $nElemT-1] } {incr i 1} {
mass [expr 4*$i+1] [expr ($Elemass($i)+$Elemass([expr $i+1]))/8] [expr ($Elemass($i)+$Elemass([expr $i+1]))/8] [expr ($Elemass($i)+$Elemass([expr $i+1]))/8]
mass [expr 4*$i+2] [expr ($Elemass($i)+$Elemass([expr $i+1]))/8] [expr ($Elemass($i)+$Elemass([expr $i+1]))/8] [expr ($Elemass($i)+$Elemass([expr $i+1]))/8]
mass [expr 4*$i+3] [expr ($Elemass($i)+$Elemass([expr $i+1]))/8] [expr ($Elemass($i)+$Elemass([expr $i+1]))/8] [expr ($Elemass($i)+$Elemass([expr $i+1]))/8]
mass [expr 4*$i+4] [expr ($Elemass($i)+$Elemass([expr $i+1]))/8] [expr ($Elemass($i)+$Elemass([expr $i+1]))/8] [expr ($Elemass($i)+$Elemass([expr $i+1]))/8]
puts $fileId " node [expr 4*$i+1]-[expr 4*$i+4] [expr ($Elemass($i)+$Elemass([expr $i+1]))/8] "
}
mass [expr 4*$nElemT+1] [expr $Elemass($nElemT)/8] [expr $Elemass($nElemT)/8] [expr $Elemass($nElemT)/8]
mass [expr 4*$nElemT+2] [expr $Elemass($nElemT)/8] [expr $Elemass($nElemT)/8] [expr $Elemass($nElemT)/8]
mass [expr 4*$nElemT+3] [expr $Elemass($nElemT)/8] [expr $Elemass($nElemT)/8] [expr $Elemass($nElemT)/8]
mass [expr 4*$nElemT+4] [expr $Elemass($nElemT)/8] [expr $Elemass($nElemT)/8] [expr $Elemass($nElemT)/8]
puts $fileId " node [expr 4*$nElemT+1]-[expr 4*$nElemT+1+4] [expr $Elemass($nElemT)/8] "
close $fileId


#-----------------------------------------------------------------------------------------
#6.CREATE SOIL ELEMENTS
#-----------------------------------------------------------------------------------------

set wgtX 0.0
set wgtZ 0.0
for {set k 1} {$k<=$numLayers} {incr k 1} {
set wgtY($k) [expr -9.81*$rho($k)]
}
set count 0
set fileId [open $dataDirSoil/SoilElements.txt w ]
# loop over layers
for {set k 1} {$k<=$numLayers} {incr k 1} {
# loop over elements
for {set j 1} {$j<=$nElemY($k)} {incr j 1} {
set nI [expr 4*($j+$count)-3]
set nJ [expr $nI+1]
set nK [expr $nI+2]
set nL [expr $nI+3]
set nM [expr $nI+4]
set nN [expr $nI+5]
set nO [expr $nI+6]
set nP [expr $nI+7]
element SSPbrick [expr $j+$count] $nI $nJ $nK $nL $nM $nN $nO $nP $k $wgtX $wgtY($k) $wgtZ
puts $fileId "[expr $j+$count] $nI $nJ $nK $nL $nM $nN $nO $nP"
}
set count [expr $count+$nElemY($k)]
}
close $fileId
puts "Finished creating all soil elements..."

#-----------------------------------------------------------------------------------------
#11. EigenValue ANALYSIS
#-----------------------------------------------------------------------------------------
# EigenValue
set numModes 20;
for { set k 1 } { $k <= $numModes } { incr k } {
recorder Node -file $dataDirSoil/mode$k-1.out -time -nodeRange 1 $nNodeT -dof 1 "eigen $k"
recorder Node -file $dataDirSoil/mode$k-2.out -time -nodeRange 1 $nNodeT -dof 1 "eigen $k"
recorder Node -file $dataDirSoil/mode$k-3.out -time -nodeRange 1 $nNodeT -dof 1 "eigen $k"
}

set lambda [eigen $numModes];

set omega {};
set Freq {};
set T {};

foreach lam $lambda {
lappend omega [expr sqrt($lam)]
lappend Freq [expr sqrt($lam)/(2*$PI)]
lappend T [expr (2*$PI)/sqrt($lam)]
}

# Outputting Freq. & Period
set Frequency [open $dataDirSoil/Freq.txt "w"]
foreach t $Freq {
puts $Frequency " $t"
}
close $Frequency

set Periods [open $dataDirSoil/period.txt "w"]
foreach t $T {
puts $Periods " $t"
}
close $Periods
puts " Analysis of EigenValue Done"
Post Reply