Problem with Domain Reduction Method for SSI

[KyriakosTryphonos]

Hello,

I am sorry for posting here, but the Soil Modeling section is locked and i can't post there.I am a post-Graduate student at Upatras, Greece and my problem is that want to model a 2D soil domain with the Domain Reduction Method (DRM). At first i want to verify the method by modeling a single Elastic 1D layer which is identical to the first example of this SFSI draft report (Nonlinear Finite Element Analysis of Dynamic Soil-Foundation-Structure Interaction) by Presisig and Jeremic 2005, which was also modeled with OpenSees.
According to the DRM method i can perform 1st analysis with the whole free field soil domain, and then record the response to a smaller boundary. Then i do not have to model all the soil domain, but only the inner boundary by i applying the boundary conditions at the inner boundary nodes. The condition that i apply in the second analysis is a ForceHistory which is found by multiplying the VelocityHistory of node with (Vs*ρsoil*Abase).The velocity History of boundary node was recorded during the whole model analysis.
Now the second analysis with the smaller domain should give identical results with the 1st one without any structure interacting with the soil, and similar results if i model a structure.
The problem is that when i try to apply this simple example, the velocity recorded at the boundary node in 2nd analysis is differs a lot from the original analysis for some reason.
Is there any possibility to provide me with a simple example of this method or give me some explanations to verify it so i can model a larger domain? Also here are the 2 scripts. I would appreciate it if you can look them and tell me if i am doing anything wrong.

Thank you for your Time.


Here is the code for the two different analysis

1)

model BasicBuilder -ndm 3 -ndf 3

set nodesInfo [open nodesInfo.dat w]

set elementInfo [open elementInfo.dat w]

set meshFile [open column.flavia.msh w]
puts $meshFile "MESH ffBrick dimension 3 ElemType Hexahedra Nnode 8"
puts $meshFile "Coordinates"

puts $meshFile "#node_number coord_x coord_y coord_z "


set LengthX 1.0
set LengthY 1.0
set LengthZ 1.0

set totLength 13.0
set nEleY 13
set nNodesY [expr $nEleY+1]






model BasicBuilder -ndm 3 -ndf 3

set CoordY 0.0

for {set i 1} {$i <= [expr $nNodesY*4]} {incr i 4} {
node $i 0.0 $CoordY 0.0
node [expr $i+1] $LengthX $CoordY 0.0
node [expr $i+2] $LengthX $CoordY $LengthZ
node [expr $i+3] 0.0 $CoordY $LengthZ



puts $nodesInfo "$i 0.0 $CoordY 0.0"
puts $nodesInfo "[expr $i+1] $LengthX $CoordY 0.0"
puts $nodesInfo "[expr $i+2] $LengthX $CoordY $LengthZ"
puts $nodesInfo "[expr $i+3] 0.0 $CoordY $LengthZ"

puts $meshFile "$i 0.0 $CoordY 0.0"
puts $meshFile "[expr $i+1] $LengthX $CoordY 0.0"
puts $meshFile "[expr $i+2] $LengthX $CoordY $LengthZ"
puts $meshFile "[expr $i+3] 0.0 $CoordY $LengthZ"

set CoordY [expr $CoordY+$LengthY]
}

for {set i 5} {$i <= [expr $nNodesY*4]} {incr i 4} {
fix $i 0 1 1
fix [expr $i+1] 0 1 1
fix [expr $i+2] 0 1 1
fix [expr $i+3] 0 1 1

}

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


for {set i 1} {$i <= [expr $nNodesY*4]} {incr i 4} {

equalDOF $i [expr $i+1] 1 2 3
equalDOF $i [expr $i+2] 1 2 3
equalDOF $i [expr $i+3] 1 2 3

}



set nNodesSoil [expr $nNodesY*4]
set c1 $nNodesSoil






puts $meshFile "end coordinates"
puts $meshFile " "
puts $meshFile "Elements"
puts $meshFile "# element node1 node2 node3 node4 node5 node6 node7 node8"


model BasicBuilder -ndm 3 -ndf 3

# poisson's ratio of soil
set PoissonSoil 0.35
# soil elastic modulus (kPa)
set Esoil 1944
# soil shear modulus (KPa)
set Gsoil [expr $Esoil/(2*(1+$PoissonSoil))]
# soil mass density (Mg/m^3)
set rho 1.8
# soil shear wave velocity (m/s)
set Vs [expr sqrt($Gsoil/$rho)]
# Element Volume
set EleVol [expr $LengthX*$LengthY*$LengthZ]
# bedrock shear wave velocity (m/s)
set rockVS 760
# bedrock mass density (Mg/m^3)
set rockDen 2.4
# Base area (m^2)
set AreaBase [expr $LengthX*$LengthZ]
# Load factor at base
set Cfactor [expr $AreaBase*$rockDen*$rockVS]

nDMaterial ElasticIsotropic 1 $Esoil $PoissonSoil $rho

set c2 1
for {set i 1} {$i <= [expr 4*$nEleY]} {incr i 4} {
element SSPbrick $c2 $i [expr $i+1] [expr $i+2] [expr $i+3] [expr $i+4] [expr $i+5] [expr $i+6] [expr $i+7] 1
puts $elementInfo "$c2 $i [expr $i+1] [expr $i+2] [expr $i+3] [expr $i+4] [expr $i+5] [expr $i+6] [expr $i+7]"
puts $meshFile "$c2 $i [expr $i+1] [expr $i+2] [expr $i+3] [expr $i+4] [expr $i+5] [expr $i+6] [expr $i+7]"
set c2 [expr $c2+1]
}

puts $meshFile "end elements"






recorder Node -file disp.out -time -nodeRange 1 [expr $c1] -dof 1 2 3 disp
recorder Node -file acceTop.out -node [expr $c1] -dof 1 accel
recorder Node -file acceBase.out -node 1 -dof 1 accel
recorder Node -file dispTop.out -time -node [expr $c1] -dof 1 disp
recorder Node -file dispBase.out -node 1 -dof 1 disp

recorder Node -file dispBound.out -node 20 -dof 1 disp
recorder Node -file velBound.out -node 20 -dof 1 vel
recorder Node -file acceBound.out -node 20 -dof 1 accel

set SineT 1.0

set pi [expr 2*asin(1)]

timeSeries Trig 1 0 5 $SineT -factor $Cfactor

pattern Plain 1 1 {
load 1 1.0 0.0 0.0
}



recorder display "OpenSees Real Time" 20 20 700 700 -wipe
prp [expr 1.0/2.0] [expr 13/2.0] 1
vup 0 1 0
vpn 0 0 1
viewWindow [expr -10] [expr 10] [expr -5] [expr 20]
display 2 10 1



integrator Newmark 0.5 0.25
numberer RCM
system FullGeneral
constraints Transformation
test NormDispIncr 1e-5 100 1
algorithm ModifiedNewton
analysis Transient

analyze 250 0.02




2)
model BasicBuilder -ndm 3 -ndf 3

set nodesInfo [open nodesInfo.dat w]

set elementInfo [open elementInfo.dat w]

set meshFile [open column.flavia.msh w]
puts $meshFile "MESH ffBrick dimension 3 ElemType Hexahedra Nnode 8"
puts $meshFile "Coordinates"

puts $meshFile "#node_number coord_x coord_y coord_z "


set LengthX 1.0
set LengthY 1.0
set LengthZ 1.0

set totLength 9.0
set nEleY 9
set nNodesY [expr $nEleY+1]






model BasicBuilder -ndm 3 -ndf 3

set CoordY 4.0

for {set i 17} {$i <= [expr $nNodesY*4+16]} {incr i 4} {
node $i 0.0 $CoordY 0.0
node [expr $i+1] $LengthX $CoordY 0.0
node [expr $i+2] $LengthX $CoordY $LengthZ
node [expr $i+3] 0.0 $CoordY $LengthZ



puts $nodesInfo "$i 0.0 $CoordY 0.0"
puts $nodesInfo "[expr $i+1] $LengthX $CoordY 0.0"
puts $nodesInfo "[expr $i+2] $LengthX $CoordY $LengthZ"
puts $nodesInfo "[expr $i+3] 0.0 $CoordY $LengthZ"

puts $meshFile "$i 0.0 $CoordY 0.0"
puts $meshFile "[expr $i+1] $LengthX $CoordY 0.0"
puts $meshFile "[expr $i+2] $LengthX $CoordY $LengthZ"
puts $meshFile "[expr $i+3] 0.0 $CoordY $LengthZ"

set CoordY [expr $CoordY+$LengthY]
}

for {set i 17} {$i <= [expr $nNodesY*4+16]} {incr i 4} {
fix $i 0 1 1
fix [expr $i+1] 0 1 1
fix [expr $i+2] 0 1 1
fix [expr $i+3] 0 1 1

}




for {set i 17} {$i <= [expr $nNodesY*4+16]} {incr i 4} {

equalDOF $i [expr $i+1] 1 2 3
equalDOF $i [expr $i+2] 1 2 3
equalDOF $i [expr $i+3] 1 2 3

}



set nNodesSoil [expr $nNodesY*4]
set c1 $nNodesSoil


# poisson's ratio of soil
set PoissonSoil 0.35
# soil elastic modulus (kPa)
set Esoil 1944
# soil shear modulus (KPa)
set Gsoil [expr $Esoil/(2*(1+$PoissonSoil))]
# soil mass density (Mg/m^3)
set rho 1.8
# soil shear wave velocity (m/s)
set Vs [expr sqrt($Gsoil/$rho)]
# Element Volume
set EleVol [expr $LengthX*$LengthY*$LengthZ]
# bedrock shear wave velocity (m/s)
set rockVS 760
# bedrock mass density (Mg/m^3)
set rockDen 2.4
# Base area (m^2)
set AreaBase [expr $LengthX*$LengthZ]
# Load factor at base
set Cfactor [expr $AreaBase*$rho*$Vs]



puts $meshFile "end coordinates"
puts $meshFile " "
puts $meshFile "Elements"
puts $meshFile "# element node1 node2 node3 node4 node5 node6 node7 node8"


model BasicBuilder -ndm 3 -ndf 3



nDMaterial ElasticIsotropic 1 $Esoil $PoissonSoil $rho

set c2 1
for {set i 17} {$i <= [expr 4*$nEleY+16]} {incr i 4} {
element SSPbrick $c2 $i [expr $i+1] [expr $i+2] [expr $i+3] [expr $i+4] [expr $i+5] [expr $i+6] [expr $i+7] 1
puts $elementInfo "$c2 $i [expr $i+1] [expr $i+2] [expr $i+3] [expr $i+4] [expr $i+5] [expr $i+6] [expr $i+7]"
puts $meshFile "$c2 $i [expr $i+1] [expr $i+2] [expr $i+3] [expr $i+4] [expr $i+5] [expr $i+6] [expr $i+7]"
set c2 [expr $c2+1]
}

puts $meshFile "end elements"






recorder Node -file disp.out -time -nodeRange 17 [expr $c1+16] -dof 1 2 3 disp
recorder Node -file acceTop1.out -node [expr $c1+16] -dof 1 accel
recorder Node -file acceBase1.out -node 20 -dof 1 accel
recorder Node -file dispTop1.out -time -node [expr $c1+16] -dof 1 disp
recorder Node -file dispBase1.out -node 20 -dof 1 disp

recorder Node -file velBase2.out -node 20 -dof 1 disp



set velHistory velBound.out
timeSeries Path 1 -dt 0.02 -filePath $velHistory -factor $Cfactor

pattern Plain 1 1 {
load 20 1.0 0.0 0.0
}



recorder display "OpenSees Real Time" 20 20 700 700 -wipe
prp [expr 1.0/2.0] [expr 13/2.0] 1
vup 0 1 0
vpn 0 0 1
viewWindow [expr -10] [expr 10] [expr -5] [expr 20]
display 2 10 1



integrator Newmark 0.5 0.25
numberer RCM
system SparseGeneral
constraints Transformation
test NormDispIncr 1e-6 400 1
algorithm KrylovNewton
analysis Transient

analyze 250 0.02


Script 2 need the velocityHistory of boundary node to work

[PeymanA]

Hey KyriakosTryphonos,

I am not sure if I understand your question correctly. In the second step, you are applying the response at the boundaries of smaller model, right?
Please help me understand your problem then we can work on it.

[KyriakosTryphonos]

Yes you are correct in the second step i am trying to apply the response as an equal force (Peff) in a smaller model with 9 elements instead of the original model which has 13 elements (thats why original Coordinate y is 4.0 in 2nd model). According to DRM theory if i apply the correct equal force in the boundary of the smaller model i should get identical response as in the 1st whole model.
You can see this small model in "Nonlinear Finite Element Analysis of Dynamic Soil-Foundation-Structure Interaction" by Presisig and Jeremic 2005 , page 16

In case of this works i can model only the smaller domain and put a structure inside to study the SSI. In my case i want to study the dynamic response of a natural gas pipe which is very long and the far field is very large.
So in step 1 i solve the far field once with no structure, and in step 2 i can model my structure with a smaller domain model with the correct boundary forces Peff.

That is the best explanation i can give. If anyone is familiar with this DRM method and can provide me some instructions i will be very grateful
Thank you!

[PeymanA]

Got it. I think your mistake is that you apply the wrong force on the boundary of small model. You should apply the response that you get from the first model minus the reflection from the surface of first model. Imagine that there is no structure, the final stress ( or force) at point i,j is the stress induced because of incident wave, minus the reflected wave after coming back from a certain surface.

You will be good if you do it this way. I encourage you to look at this thesis available at OpenSees website.
http://opensees.berkeley.edu/wiki/index ... ad_Pattern


Hope that helps.

[KyriakosTryphonos]

Thank you very much.
I understand now what i was doing wrong and i will try to fix it. You 've been very helpful!