Eigen value for saturated soil column

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

Moderators: silvia, selimgunay, Moderators

Post Reply
masoumeh
Posts: 5
Joined: Thu Feb 18, 2016 11:32 pm

Eigen value for saturated soil column

Post by masoumeh » Tue Nov 22, 2022 6:42 pm

Hi,

I want to find the natural period of the saturated soil column with coupled soil-fluid elements. When I am modeling dry soil column, I can get the eigen values, but when I add pore water pressure dof and use coupled soil-fluid elements, use of command "eigen 1" gives me error. The resulting stress and strains of my saturated coupled model is correct, I only cannot do eigen analysis. Can anyone help me how I can address this problem?

I am using 3D model with 4 dofs. I tried different analyze parameters, but the problem is not resolved.

constraints Penalty 1.0e14 1.0e14
test NormDispIncr 1.0e-3 1000 0
integrator HHT 0.67
algorithm KrylovNewton
system SparseGEN
numberer Plain
rayleigh $a0 0.0 $a1 0.0
analysis Transient

Error: WARNING SuperLU::solve(void)- Error 4 returned in factorization dgstrf()

masoumeh
Posts: 5
Joined: Thu Feb 18, 2016 11:32 pm

Re: Eigen value for saturated soil column

Post by masoumeh » Tue Nov 22, 2022 6:46 pm

Here is the code:

set alfa 4.0
# size of model, unit is meter (m)
set modelX 0.5
set modelZ 10.0
set modelY 0.5
# number of elements
set nElemX 1
set nElemZ 20
set nElemY 1
# water table being relative to surface, unit is meter (m)
set wt 1.0

# Parameters of materials and thickness of each layers
set Dr 0.55
set Dr_Nliq 0.9
set H1 1.5
set H2 4.0
set H3 [expr 10-$H1-$H2]

############################################################
# size of element (on the center axis)
set sElemX [expr $modelX/$nElemX]
set sElemY [expr $modelY/$nElemY]
set sElemZ [expr $modelZ/$nElemZ]

set nNodeX [expr $nElemX+1]
set nNodeY [expr $nElemY+1]
set nNodeZ [expr $nElemZ+1]

set nElemT [expr $nElemX*$nElemY*$nElemZ]
set nNodeL [expr $nNodeX*$nNodeY]
set nNodeT [expr $nNodeZ*$nNodeL]

set Gini 2.0e5
set poisini 0.35
set g [expr -1.0 * $ng]

set dur [expr $motionSteps*$dms/1000.0]

#-----------------------------------------------------------------------------------------
# 2. CREATE SOIL NODES AND BOUNDARY CONDITIONS
#-----------------------------------------------------------------------------------------

model BasicBuilder -ndm 3 -ndf 4

set nodesInfo [open nodesInfo.dat w]

for {set i 1} {$i <= $nNodeZ} {incr i} {
for {set j 1} {$j <= $nNodeY} {incr j} {
for {set k 1} {$k <= $nNodeX} {incr k} {
set posX [expr ($k-1)*$sElemX]
set posY [expr ($j-1)*$sElemY]
set posZ [expr ($i-1)*$sElemZ]
set id [expr ($i-1)*$nNodeL+($j-1)*$nNodeX+$k]
# puts "$id $posX $posY $posZ"
node $id $posX $posY $posZ
puts $nodesInfo "$id $posX $posY $posZ"

#mass $id 0.00000000000001 0.00000000000001 0.00000000000001 0.00000000000001
set f1 0
set f2 1
set f3 0
set f4 0
if {$id == 1} {
set f1 1
set f3 1
}
if {$posZ >= $modelZ-$wt} {
set f4 1
}
fix $id $f1 $f2 $f3 $f4

set equalid [expr ($i-1)*$nNodeL+1]
if {$id > $equalid} {
equalDOF $equalid $id 1 3
}
}
}
}

close $nodesInfo
puts "Finished creating all soil nodes..."


#-----------------------------------------------------------------------------------------
# 4. CREATE SOIL MATERIALS
#-----------------------------------------------------------------------------------------

# from bottom layer to top layer
# the "ytopn" means distance from the bottom of soil column to the top of one layer
set e_max 0.78
set e_min 0.51
set e_in [expr $e_max-$Dr*($e_max-$e_min)]
puts "e_in= $e_in"
set e_in_Nliq [expr $e_max-$Dr_Nliq*($e_max-$e_min)]
puts "e_in_Nliq= $e_in_Nliq"
# bottom layer
set eini1 $e_in_Nliq
set ytop1 $H3
# middle layer
set eini2 $e_in
set ytop2 [expr $H3+$H2]
# top layer
set eini3 $e_in_Nliq
set ytop3 10.0
###########################

set den1 [expr 2.65/($eini1+1.0)*(1+$eini1/2.65)]
set den2 [expr 2.65/($eini2+1.0)*(1+$eini2/2.65)]
set den3 [expr 2.65/($eini3+1.0)*(1+$eini3/2.65)]
set den4 [expr 2.65/($eini3+1.0)]
puts "den1 = $den1"
puts "den2 = $den2"
puts "den3 = $den3"
puts "den4 = $den4"

# initial void ratio
set pois 0.05

# SanisandMSf tag G0 nu eini Mc c lambda_c e0
# ksi P_atm m h0pie ch nsupb A0pie nsupd
# z_max cz nsubg miu0 u cl x cr
# Den integMode TolR psh pcut model
# psh = threshold for the semi-fluidized state
# model: 0 = elastic; 1 = Sanisand-MSf; 2 = MD04
# integMode: 1 = modified Euler with error control; 2 = forward; 3 = backward

nDMaterial SanisandMSf2 1 125 $pois $eini1 1.26 0.735 0.0287 0.78 \
0.7 101 0.01 4.6 0.968 3.5 0.626 2.5 \
15 2000 0.9 1.99 1.32 35 3.5 0.0 \
$den1 1 1e-4 10 0.5 1
nDMaterial SanisandMSf2 2 125 $pois $eini2 1.26 0.735 0.0287 0.78 \
0.7 101 0.01 4.6 0.968 3.5 0.626 2.5 \
15 2000 0.9 1.99 1.32 35 3.5 0.0 \
$den2 1 1e-4 10 0.5 1
nDMaterial SanisandMSf2 3 125 $pois $eini3 1.26 0.735 0.0287 0.78 \
0.7 101 0.01 4.6 0.968 3.5 0.626 2.5 \
15 2000 0.9 1.99 1.32 35 3.5 0.0 \
$den3 1 1e-4 10 0.5 1
nDMaterial SanisandMSf2 4 125 $pois $eini1 1.26 0.735 0.0287 0.78 \
0.7 101 0.01 4.6 0.968 3.5 0.626 2.5 \
15 2000 0.9 1.99 1.32 35 3.5 0.0 \
$den4 1 1e-4 10 0.5 1

puts "Finished creating all soil materials..."

#-----------------------------------------------------------------------------------------
# 5. CREATE SOIL ELEMENTS
#-----------------------------------------------------------------------------------------

set elemInfo [open elementInfo.dat w]
# loop over layers

set facx [expr sin($alfa*3.1415/180)]
set facy [expr -cos($alfa*3.1415/180)]
set gx [expr $ng*$facx]
set gy [expr $ng*$facy]
set uBulk 2.0e6
set alpha 1e-7

for {set i 1} {$i <= $nElemZ} {incr i} {
for {set j 1} {$j <= $nElemY} {incr j} {
for {set k 1} {$k <= $nElemX} {incr k} {
set nI [expr ($i-1)*$nNodeL+($j-1)*$nNodeX+$k]
set nJ [expr ($i-1)*$nNodeL+($j-1)*$nNodeX+$k+1]
set nK [expr ($i-1)*$nNodeL+$j*$nNodeX+$k+1]
set nL [expr ($i-1)*$nNodeL+$j*$nNodeX+$k]
set nM [expr $i*$nNodeL+($j-1)*$nNodeX+$k]
set nN [expr $i*$nNodeL+($j-1)*$nNodeX+$k+1]
set nO [expr $i*$nNodeL+$j*$nNodeX+$k+1]
set nP [expr $i*$nNodeL+$j*$nNodeX+$k]

set id [expr ($i-1)*$nElemX*$nElemY+($j-1)*$nElemX+$k]
# puts "$id $nI $nJ $nK $nL $nM $nN $nO $nP"
set posZ [expr ($i-0.5)*$sElemZ]
set matNum 1
set Den $den1
set eini $eini1

if {$posZ <= $ytop1} {
set matNum 1
set Den $den1
set eini $eini1
} elseif {$posZ <= $ytop2} {
set matNum 2
set Den $den2
set eini $eini2
} elseif {$posZ <= $modelZ-$wt} {
set matNum 3
set Den $den3
set eini $eini3
} else {
set matNum 4
set Den $den4
set eini $eini3
}
set perm [expr (0.0207*$eini-0.0009)/9.81*0.01]
element SSPbrickUP $id $nI $nJ $nK $nL $nM $nN $nO $nP $matNum \
$uBulk 1.0 $perm $perm $perm $eini $alpha $gx 0.0 $gy
#element SSPbrick $id $nI $nJ $nK $nL $nM $nN $nO $nP $matNum $gx 0.0 $gy

puts $elemInfo "$id $nI $nJ $nK $nL $nM $nN $nO $nP"
}
}
}


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

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

set conNum 10
set conDT [expr $conNum*$dT]

recorder Element -file output/Gstress.out -time -dT $conDT -eleRange 1 $nElemT stress
recorder Node -file output/Gdisp.out -time -dT $conDT -nodeRange 1 $nNodeT -dof 1 2 3 disp
recorder Node -file output/Gpwp.out -time -dT $conDT -nodeRange 1 $nNodeT -dof 4 vel

puts "Finished creating gravity recorders..."


#-----------------------------------------------------------------------------------------
# 7. GRAVITY ANALYSIS
#-----------------------------------------------------------------------------------------


model BasicBuilder -ndm 3 -ndf 3

# Rayleigh damping parameter
# in flac3d, the damp ratio=1%, center frequency=1Hz
set fmin 2.5
set dampRatio 0.01
set a1 [expr $dampRatio/2.0/3.1415/$fmin]
set a0 [expr $dampRatio*2.0*3.1415*$fmin]

logFile Log.txt
# Create analysis
constraints Penalty 1.0e14 1.0e14
test NormDispIncr 1.0e-3 1000 0
#integrator Newmark 0.5 0.25
integrator HHT 0.67
algorithm KrylovNewton ;#
system SparseGEN ;# Use sparse solver. Next numberer is better to be Plain.
numberer Plain ;# method to map between between equation numbers of DOFs

rayleigh $a0 0.0 $a1 0.0
analysis Transient

set startT [clock seconds]


# gravity loading
puts "elastic consolidation start..."
setParameter -value $poisini -eleRange 1 $nElemT poissonRatio 1
setParameter -value $poisini -eleRange 1 $nElemT poissonRatio 2
setParameter -value $poisini -eleRange 1 $nElemT poissonRatio 3
setParameter -value $poisini -eleRange 1 $nElemT poissonRatio 4
setParameter -value $Gini -eleRange 1 $nElemT fixGforS0 1
setParameter -value $Gini -eleRange 1 $nElemT fixGforS0 2
setParameter -value $Gini -eleRange 1 $nElemT fixGforS0 3
setParameter -value $Gini -eleRange 1 $nElemT fixGforS0 4
updateMaterialStage -material 1 -stage 0
updateMaterialStage -material 2 -stage 0
updateMaterialStage -material 3 -stage 0
updateMaterialStage -material 4 -stage 0

analyze 100 1

setParameter -value $pois -eleRange 1 $nElemT poissonRatio 1
setParameter -value $pois -eleRange 1 $nElemT poissonRatio 2
setParameter -value $pois -eleRange 1 $nElemT poissonRatio 3
setParameter -value $pois -eleRange 1 $nElemT poissonRatio 4
analyze 100 1

puts "calm down"
updateMaterialStage -material 1 -stage 1
updateMaterialStage -material 2 -stage 1
updateMaterialStage -material 3 -stage 1
updateMaterialStage -material 4 -stage 1

# analyze 100 $dT

puts "Finished creating all soil boundary conditions..."
remove recorders


#set eigenvalues [eigen -fullGenLapack 10];
set eigenvalues [eigen -genBandArpack 5];
puts $eigenvalues
wipe

Post Reply