problems with u-p-U & u-p elements

[Shahir]

Dear Moderator and users
I have some problems with coupled elements (u-p-U & u-p elements) as follows:

1. I analyzed 1D consolidation and triaxial model examples using u-p-U element (the input files have been attached). The results of version 1.6.0 for consolidation example are wrong but the version 1.7.2 produces correct results. Vice versa, In the case of triaxial example, the results are correct using version 1.6.0 and wrong using version 1.7.2. I really don’t know in which version u-p-U element works correctly?

2. I analyzed 1D consolidation using u-p element. Both of versions 1.6.0 & 1.7.2 produce unfavorable results. Does this element work?

Also, I have a question:

3. How Impermeable boundary can be modeled using u-p-U model. You know that the relative velocity of fluid phase shall be fixed in the impermeable boundary (or in the other words u = U). But, there is no option to fix relative velocity of fluid phase.

Please, let me to know your idea. Thank you in forward for your help.
Sincerely
Shahir



#========================#
# 1D Consolidation _ upU element #
#========================#

wipe all

#Basic units
# meter for length
set m 1.0
# second for time
set sec 1.0
# Kilogram for mass
set kg 1.0

# angle
set rad 1.0
set PI [expr 4*atan(1.0)]
set deg [expr $PI/180.0*$rad]

# length
set cm [expr 0.01*$m]
set in 0.0254
set ft [expr 12.0*$in]

# mass
set lbs 0.4536
set kip [expr 1000.0*$lbs]
set ton [expr 1000.0*$kg]

# force
set N [expr 1.0*$kg*$m/pow($sec,2)]
set kN [expr 1000.0*$N]

# pressure
set Pa [expr 1.0*$N/pow($m, 2)]
set kPa [expr 1000.0*$Pa]
set pcf [expr $lbs/pow($ft,3)]
set ksi [expr $kip/pow($in,2)]
set psi [expr $ksi/1000.]

# Create the modelbuilder
model BasicBuilder -ndm 3 -ndf 7

# Nodal coordinates
# tag X Y Z
node 1 0.5 0.5 0.0
node 2 0.0 0.5 0.0
node 3 0.0 0.0 0.0
node 4 0.5 0.0 0.0
node 5 0.5 0.5 0.5
node 6 0.0 0.5 0.5
node 7 0.0 0.0 0.5
node 8 0.5 0.0 0.5
node 9 0.5 0.5 1.0
node 10 0.0 0.5 1.0
node 11 0.0 0.0 1.0
node 12 0.5 0.0 1.0
node 13 0.5 0.5 1.5
node 14 0.0 0.5 1.5
node 15 0.0 0.0 1.5
node 16 0.5 0.0 1.5
node 17 0.5 0.5 2.0
node 18 0.0 0.5 2.0
node 19 0.0 0.0 2.0
node 20 0.5 0.0 2.0
node 21 0.5 0.5 2.5
node 22 0.0 0.5 2.5
node 23 0.0 0.0 2.5
node 24 0.5 0.0 2.5
node 25 0.5 0.5 3.0
node 26 0.0 0.5 3.0
node 27 0.0 0.0 3.0
node 28 0.5 0.0 3.0
node 29 0.5 0.5 3.5
node 30 0.0 0.5 3.5
node 31 0.0 0.0 3.5
node 32 0.5 0.0 3.5
node 33 0.5 0.5 4.0
node 34 0.0 0.5 4.0
node 35 0.0 0.0 4.0
node 36 0.5 0.0 4.0

# Boundary conditions
# nodeID ux uy uz p Ux Uy Uz
fix 1 1 1 1 0 1 1 1
fix 2 1 1 1 0 1 1 1
fix 3 1 1 1 0 1 1 1
fix 4 1 1 1 0 1 1 1
fix 5 1 1 0 0 1 1 0
fix 6 1 1 0 0 1 1 0
fix 7 1 1 0 0 1 1 0
fix 8 1 1 0 0 1 1 0
fix 9 1 1 0 0 1 1 0
fix 10 1 1 0 0 1 1 0
fix 11 1 1 0 0 1 1 0
fix 12 1 1 0 0 1 1 0
fix 13 1 1 0 0 1 1 0
fix 14 1 1 0 0 1 1 0
fix 15 1 1 0 0 1 1 0
fix 16 1 1 0 0 1 1 0
fix 17 1 1 0 0 1 1 0
fix 18 1 1 0 0 1 1 0
fix 19 1 1 0 0 1 1 0
fix 20 1 1 0 0 1 1 0
fix 21 1 1 0 0 1 1 0
fix 22 1 1 0 0 1 1 0
fix 23 1 1 0 0 1 1 0
fix 24 1 1 0 0 1 1 0
fix 25 1 1 0 0 1 1 0
fix 26 1 1 0 0 1 1 0
fix 27 1 1 0 0 1 1 0
fix 28 1 1 0 0 1 1 0
fix 29 1 1 0 0 1 1 0
fix 30 1 1 0 0 1 1 0
fix 31 1 1 0 0 1 1 0
fix 32 1 1 0 0 1 1 0
fix 33 1 1 0 1 1 1 0
fix 34 1 1 0 1 1 1 0
fix 35 1 1 0 1 1 1 0
fix 36 1 1 0 1 1 1 0

# Parameters
set grvt [expr 9.8*$m/pow($sec,2)]
set E [expr 1.0e4*$kPa]
set v 0.35
set poro 0.4
set alpha 1.0
set rho_s [expr 2.0e3*$kg/pow($m,3)]
set rho_f [expr 1.0e3*$kg/pow($m,3)]
set bulk_f [expr 2.0e9*$Pa]
set bulk_s [expr 3.6e10*$Pa]
set kx [expr 1.0e-4*$m/$sec/($grvt*$rho_f)]
set ky [expr 1.0e-4*$m/$sec/($grvt*$rho_f)]
set kz [expr 1.0e-4*$m/$sec/($grvt*$rho_f)]

# Material type and Elements
# tag E nu density
nDMaterial ElasticIsotropic3D 1 $E $v $rho_s

# tag <-------8 Nodes-------> matID BF1 BF2 BF3 porosity alpha solid_density fluid_density perm_x perm_y perm_z s_bulk_modu f_bulk_modu pressure
element Brick8N_u_p_U 1 1 2 3 4 5 6 7 8 1 0.0 0.0 0.0 $poro $alpha $rho_s $rho_f $kx $ky $kz $bulk_s $bulk_f 0.0
element Brick8N_u_p_U 2 5 6 7 8 9 10 11 12 1 0.0 0.0 0.0 $poro $alpha $rho_s $rho_f $kx $ky $kz $bulk_s $bulk_f 0.0
element Brick8N_u_p_U 3 9 10 11 12 13 14 15 16 1 0.0 0.0 0.0 $poro $alpha $rho_s $rho_f $kx $ky $kz $bulk_s $bulk_f 0.0
element Brick8N_u_p_U 4 13 14 15 16 17 18 19 20 1 0.0 0.0 0.0 $poro $alpha $rho_s $rho_f $kx $ky $kz $bulk_s $bulk_f 0.0
element Brick8N_u_p_U 5 17 18 19 20 21 22 23 24 1 0.0 0.0 0.0 $poro $alpha $rho_s $rho_f $kx $ky $kz $bulk_s $bulk_f 0.0
element Brick8N_u_p_U 6 21 22 23 24 25 26 27 28 1 0.0 0.0 0.0 $poro $alpha $rho_s $rho_f $kx $ky $kz $bulk_s $bulk_f 0.0
element Brick8N_u_p_U 7 25 26 27 28 29 30 31 32 1 0.0 0.0 0.0 $poro $alpha $rho_s $rho_f $kx $ky $kz $bulk_s $bulk_f 0.0
element Brick8N_u_p_U 8 29 30 31 32 33 34 35 36 1 0.0 0.0 0.0 $poro $alpha $rho_s $rho_f $kx $ky $kz $bulk_s $bulk_f 0.0

# Loading
set p [expr 100*$kN]
set sin "Sine 0 1e-5 4e-5"
set rec "Rectangular 1e-5 1000"
pattern Plain 1 $sin {
# nodeID Fx Fy Fz
load 33 0 0 $p 0 0 0 0
load 34 0 0 $p 0 0 0 0
load 35 0 0 $p 0 0 0 0
load 36 0 0 $p 0 0 0 0
}
pattern Plain 2 $rec {
# nodeID Fx Fy Fz
load 33 0 0 $p 0 0 0 0
load 34 0 0 $p 0 0 0 0
load 35 0 0 $p 0 0 0 0
load 36 0 0 $p 0 0 0 0
}

# Recorder generation
recorder Node -file node4 -time -node 4 -dof 3 4 disp
recorder Node -file node12 -time -node 12 -dof 3 4 disp
recorder Node -file node20 -time -node 20 -dof 3 4 disp
recorder Node -file node28 -time -node 28 -dof 3 4 disp
recorder Node -file node36 -time -node 36 -dof 3 4 disp

# Create the analysis
integrator Newmark 0.6 0.3025
test NormDispIncr 1e-6 20 1
algorithm Newton
constraints Penalty 1e12 1e12
numberer RCM
system UmfPack
analysis Transient

# Perform the analysis
puts \n
puts "Analysis begins"
set dt 1.0
puts \n
for {set i 1} {$i <=130} {incr i} {
puts -nonewline stdout "Runing step $i \b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b"
flush stdout
analyze 1 $dt
}
puts \n
puts "Analysis ends"
puts \n

wipe




#==============================#
# Triaxial model _ upU element #
#==============================#

wipe all

#Basic units
# meter for length
set m 1.0
# second for time
set sec 1.0
# Kilogram for mass
set kg 1.0

# angle
set rad 1.0
set PI [expr 4*atan(1.0)]
set deg [expr $PI/180.0*$rad]

# length
set cm [expr 0.01*$m]
set in 0.0254
set ft [expr 12.0*$in]

# mass
set lbs 0.4536
set kip [expr 1000.0*$lbs]
set ton [expr 1000.0*$kg]

# force
set N [expr 1.0*$kg*$m/pow($sec,2)]
set kN [expr 1000.0*$N]

# pressure
set Pa [expr 1.0*$N/pow($m, 2)]
set kPa [expr 1000.0*$Pa]
set pcf [expr $lbs/pow($ft,3)]
set ksi [expr $kip/pow($in,2)]
set psi [expr $ksi/1000.]
#Parameters
set grvt [expr 9.8*$m/pow($sec,2)]
set E [expr 1.2e7*$Pa]
set v 0.35
set poro 0.4
set alpha 1.0
set rho_s [expr 2.7e3*$kg/pow($m,3)]
set rho_f [expr 1.0e3*$kg/pow($m,3)]
set bulk_f [expr 2.2e9*$Pa]
set bulk_s [expr 1.0e15*$Pa]
set kx [expr 1.0e-10*$m/$sec/($grvt*$rho_f)]
set ky [expr 1.0e-10*$m/$sec/($grvt*$rho_f)]
set kz [expr 1.0e-10*$m/$sec/($grvt*$rho_f)]

#Create the modelbuilder
model BasicBuilder -ndm 3 -ndf 7

#Nodal coordinates
node 1 [expr 1.0*$m] [expr 1.0*$m] [expr 0.0*$m]
node 2 [expr 0.0*$m] [expr 1.0*$m] [expr 0.0*$m]
node 3 [expr 0.0*$m] [expr 0.0*$m] [expr 0.0*$m]
node 4 [expr 1.0*$m] [expr 0.0*$m] [expr 0.0*$m]
node 5 [expr 1.0*$m] [expr 1.0*$m] [expr 1.0*$m]
node 6 [expr 0.0*$m] [expr 1.0*$m] [expr 1.0*$m]
node 7 [expr 0.0*$m] [expr 0.0*$m] [expr 1.0*$m]
node 8 [expr 1.0*$m] [expr 0.0*$m] [expr 1.0*$m]

#Boundary conditions
fix 1 0 0 1 1 0 0 0
fix 2 0 0 1 1 0 0 0
fix 3 0 0 1 1 0 0 0
fix 4 0 0 1 1 0 0 0
fix 5 0 0 0 1 0 0 0
fix 6 0 0 0 1 0 0 0
fix 7 0 0 0 1 0 0 0
fix 8 0 0 0 1 0 0 0

# elastic material
nDMaterial ElasticIsotropic3D 1 $E $v $rho_s

# Drucker-Prager material
set FA 20 ; #Friction Angle
set DA 20 ; #Dilation Angle
set FAr [expr $FA*$PI/180]
set DAr [expr $DA*$PI/180]
set a1 [expr 2*sin($FAr)/sqrt(3)/(3-sin($FAr))]
set a2 [expr 2*sin($DAr)/sqrt(3)/(3-sin($DAr))]
set k 0.0
set ys "-DP" ; # Yield function
set ps "-DP $a2" ; # Potencial function (Flow rule)
set es1 "-Leq 0.0" ; # Isotropic hardening
set es2 "-Leq 0.0" ; # Isotropic hardening
set et1 "-Linear 0.0" ; # Kinematic hardening
set strs "-0.01 0 0 0 -0.01 0 0 0 -0.01" ; #Initial Stress
set eps "-NOD 0 -NOS 2 $a1 $k -stressp $strs"
nDMaterial Template3Dep 2 1 -YS $ys -PS $ps -EPS $eps

#_____________________tag_____8 nodes_____matID__bforce1_bforce2_bforce3 porosity alpha solid_density fluid_density perm_x perm_y perm_z s_bulk_modulus f_bulk_modulus pressure
element Brick8N_u_p_U 1 3 4 1 2 7 8 5 6 2 0.0 0.0 0.0 $poro $alpha $rho_s $rho_f $kx $ky $kz $bulk_s $bulk_f 0.0

#Isotropic Loading#
set mp -25
set np 25
set rect "Rectangular 0 1000"
pattern Plain 1 $rect {
load 1 [expr $np*$kPa] [expr $np*$kPa] [expr $mp*$kPa] 0 0 0 0
load 2 [expr $mp*$kPa] [expr $np*$kPa] [expr $mp*$kPa] 0 0 0 0
load 3 [expr $mp*$kPa] [expr $mp*$kPa] [expr $mp*$kPa] 0 0 0 0
load 4 [expr $np*$kPa] [expr $mp*$kPa] [expr $mp*$kPa] 0 0 0 0
load 5 [expr $np*$kPa] [expr $np*$kPa] [expr $np*$kPa] 0 0 0 0
load 6 [expr $mp*$kPa] [expr $np*$kPa] [expr $np*$kPa] 0 0 0 0
load 7 [expr $mp*$kPa] [expr $mp*$kPa] [expr $np*$kPa] 0 0 0 0
load 8 [expr $np*$kPa] [expr $mp*$kPa] [expr $np*$kPa] 0 0 0 0}

recorder Element -file stress_a -time -ele 1 stresses
recorder Node -file node1_a -time -node 1 -dof 1 2 3 4 5 6 7 disp
recorder Node -file node5_a -time -node 5 -dof 1 2 3 4 5 6 7 disp

integrator Newmark 0.6 0.3025
numberer RCM
constraints Penalty 1e12 1e12
test NormDispIncr 1.0e-6 50 1
algorithm Newton
system UmfPack
analysis Transient

analyze 1 1
#Axial Loading

loadConst -time 5

set d [expr -5.0e-4*$m]
pattern Plain 2 Linear {
sp 5 3 $d
sp 6 3 $d
sp 7 3 $d
sp 8 3 $d}

recorder Element -file stress_b -time -ele 1 stresses
recorder Node -file node1_b -time -node 1 -dof 1 2 3 4 5 6 7 disp
recorder Node -file node5_b -time -node 5 -dof 1 2 3 4 5 6 7 disp

analyze 50 1




#=======================#
# 1D Consolidation _ up element #
#=======================#

wipe all

#Basic units
# meter for length
set m 1.0
# second for time
set sec 1.0
# Kilogram for mass
set kg 1.0

# angle
set rad 1.0
set PI [expr 4*atan(1.0)]
set deg [expr $PI/180.0*$rad]

# length
set cm [expr 0.01*$m]
set in 0.0254
set ft [expr 12.0*$in]

# mass
set lbs 0.4536
set kip [expr 1000.0*$lbs]
set ton [expr 1000.0*$kg]

# force
set N [expr 1.0*$kg*$m/pow($sec,2)]
set kN [expr 1000.0*$N]

# pressure
set Pa [expr 1.0*$N/pow($m, 2)]
set kPa [expr 1000.0*$Pa]

# Create the modelbuilder
model BasicBuilder -ndm 2 -ndf 3

# Nodal coordinates
for {set i 1} {$i <=17} {incr i 2} {
# tag X Y
node $i 0.0 [expr ($i-1)*0.25]
}
for {set i 2} {$i <=18} {incr i 2} {
# tag X Y
node $i 0.5 [expr ($i-2)*0.25]
}

# Boundary conditions
# nodeID ux uy p
fix 1 1 1 0
fix 2 1 1 0
fix 3 1 0 0
fix 4 1 0 0
fix 5 1 0 0
fix 6 1 0 0
fix 7 1 0 0
fix 8 1 0 0
fix 9 1 0 0
fix 10 1 0 0
fix 11 1 0 0
fix 12 1 0 0
fix 13 1 0 0
fix 14 1 0 0
fix 15 1 0 0
fix 16 1 0 0
fix 17 1 0 1
fix 18 1 0 1

# Parameters
set grvt [expr 9.8*$m/pow($sec,2)]
set E [expr 1.0e4*$kPa]
set v 0.35
set rho_s [expr 2.0e3*$kg/pow($m,3)]
set rho_f [expr 1.0e3*$kg/pow($m,3)]
set bulk [expr 5.0e9*$Pa]
set kx [expr 1.0e-4*$m/$sec/($grvt*$rho_f)]
set ky [expr 1.0e-4*$m/$sec/($grvt*$rho_f)]

# Material type and Elements
# tag E nu density
nDMaterial ElasticIsotropic3D 1 $E $v $rho_s

# tag <-4 Nodes-> thick matID bulk_modu fluid_density perm_h perm_v BF1 BF2 traction
element quadUP 1 1 2 4 3 1 1 $bulk $rho_f $kx $ky 0.0 0.0 0.0
element quadUP 2 3 4 6 5 1 1 $bulk $rho_f $kx $ky 0.0 0.0 0.0
element quadUP 3 5 6 8 7 1 1 $bulk $rho_f $kx $ky 0.0 0.0 0.0
element quadUP 4 7 8 10 9 1 1 $bulk $rho_f $kx $ky 0.0 0.0 0.0
element quadUP 5 9 10 12 11 1 1 $bulk $rho_f $kx $ky 0.0 0.0 0.0
element quadUP 6 11 12 14 13 1 1 $bulk $rho_f $kx $ky 0.0 0.0 0.0
element quadUP 7 13 14 16 15 1 1 $bulk $rho_f $kx $ky 0.0 0.0 0.0
element quadUP 8 15 16 18 17 1 1 $bulk $rho_f $kx $ky 0.0 0.0 0.0

# Loading
set p [expr -400*$kN]
set sin "Sine 0 1e-5 4e-5"
set rec "Rectangular 1e-5 10000"
pattern Plain 1 $sin {
# nodeID Fx Fy
load 17 0 $p 0
load 18 0 $p 0
}
pattern Plain 2 $rec {
# nodeID Fx Fy
load 17 0 $p 0
load 18 0 $p 0
}

# Recorder generation
recorder Node -file 1_node2 -time -node 2 -dof 2 3 disp
recorder Node -file 1_node6 -time -node 6 -dof 2 3 disp
recorder Node -file 1_node10 -time -node 10 -dof 2 3 disp
recorder Node -file 1_node14 -time -node 14 -dof 2 3 disp
recorder Node -file 1_node18 -time -node 18 -dof 2 3 disp

# Create the analysis
integrator Newmark 0.6 0.3025
test NormDispIncr 1e-6 20 1
algorithm Newton
constraints Penalty 1e12 1e12
numberer RCM
system UmfPack
analysis Transient

# Perform the analysis
puts \n
puts "Analysis begins"
set dt 0.1
puts \n
for {set i 1} {$i <=1300} {incr i} {
puts -nonewline stdout "Runing step $i \b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b"
flush stdout
analyze 1 $dt
}
puts \n
puts "Analysis ends"
puts \n

wipe

[Boris]

Hello Shahir,

I am not sure about versioning of OpenSees, do not have contrlo over that (can you give me dates of versions), but we do have examples for consolidation (both small and large deformation) that works quite well. I might show some of those at the workshop in few weeks.

WE quite a bit of new development which has not made it to the OpenSees repository yet (but is available from our subversion system, see my web site)...

We also do have our own implementation of u-p (the one you used is most likely from UCSD) and this u-p (ours) also works well for consolidation...

As for boundary conditions for u-p-U, we do have a command that will link u and U DOFs, and can be used for impermeable conditions. It will make it to OpenSees CVS soon (hopefully before the workshop).

Boris

[Shahir]

Dear Professor Jeremic
Thanks for your attention.
I am working with the last version of OpenSees (1.7.2). I analyzed some examples with upU element and elastic material and there was no problem. My problem is about drained triaxial model with Drucker-Prager material. In this model, when the stress state reaches the yield surface, the confining pressure increases and the element remains in the elastic state. However, using Brick element instead of upU element, the results of this model become favorable.
It is nice of you to help me.
Best regards
Shahir

[zcheng]

Hello Shahir,

Some updates on upU element are uploaded today. Try this new version. Also, I suggest to use NewTemplate3Dep to specify the Drucker-Prager model.

Regards,
Z. Cheng

[Shahir]

Dear Cheng
Thanks for your information.
Best regards
Shahir

[macarojaime]

Professor Jeremic,

You mentioned some time ago (in your reply to Shahir on Wed Aug 02, 2006) that you have a command to link the DOFs of u and U in your upU elements. Could you please provide me with a reference to find this command in Opensees?

I am working on modeling Verdugo and Ishihara's undrained triaxial tests, and one of the essential features is to model solid and fluid parts moving simultaneously in response to an applied load.

On the other hand, I've been trying to figure out the way of keeping a constant lateral pressure throughout the shearing stage of the CU triaxial, without success. Could you recommend me something in this regards?

I would appreciate any comment.

Thank you,

MARIA JAIME
Doctoral Student
Dept. of Civil and Environm. Eng.
University of Pittsburgh