PressureDependMultiYield-quadUP element

Command_Manual
- Tcl Commands
- Modeling_Commands
  - model
  - uniaxialMaterial
  - ndMaterial
  - frictionModel
  - section
  - geometricTransf
  - element
  - node
  - sp commands
    - fix
    - fixX
    - fixY
    - fixZ
  - mp commands
  - timeSeries
  - pattern
  - mass
  - block commands
    - block2D
    - block3D
  - region
  - rayleigh
- Analysis Commands
- Output Commands
  - recorder
  - record
  - print
  - printA
  - logFile
  - RealTime_Output_Commands
- Misc Commands
- DataBase Commands
  - Template:DataBaseCommandsMenu

Solid-fluid fully coupled (u-p) plane-strain element: Inclined (4 degrees), saturated soil column with pressure dependent material subjected to 1D sinusoidal base shaking

Input File

#Created by Zhaohui Yang (zhyang@ucsd.edu)
#
#fully coupled, u-p formulation
#
#plane strain,  shear-beam type mesh with single material,  
#dynamic analysis,  SI units (m, s, KN, ton)
#
wipe
#
#some user defined variables
# 
set matOpt   1      ;# 1 = pressure depend; 
                    ;# 2 = pressure independ; 
set fmass  1      ;# fluid mass density
set smass  2      ;# saturated soil mass density
set G     6.e4     ;
set B     2.4e5    ;
set bulk   2.2e6  ;#fluid-solid combined bulk modulus
set vperm  5.e-4  ;#vertical permeability (m/s)
set hperm  5.e-4  ;#horizontal permeability (m/s)
set accGravity  9.81    ;#acceleration of gravity
set vperm  [expr $vperm/$accGravity/$fmass]  ;#actual value used in computation
set hperm  [expr $hperm/$accGravity/$fmass]  ;#actual value used in computation
set press   0.    ;# isotropic consolidation pressure on quad element(s)
set loadBias 0.07     ;# Static shear load, in percentage of gravity load (=sin(inclination angle))

set accMul  2.  ;# acc. multiplier
set accNam  whatever.acc  ;# file name for input acc. record 
set accDt   0.0166   ;# dt of input acc. record
set period   1.0      ;# Period for applied Sine wave
set deltaT   0.01    ;# time step for analysis
set numSteps 2400    ;# number of time steps
set gamma    0.6    ;# Newmark integration parameter

set massProportionalDamping   0. ;
set InitStiffnessProportionalDamping 0.002;

set numXele 1      ;# number of elements in x (H) direction
set numYele 10      ;# number of elements in y (V) direction
set xSize 1.0       ;# x direction element size
set ySize 1.0       ;# y direction element size

#############################################################
# BUILD MODEL

#create the ModelBuilder
model basic -ndm 2 -ndf 3

# define material and properties
switch $matOpt {
  1 {
    nDMaterial PressureDependMultiYield 1 2 $smass $G $B 31.4 .1 80 0.5\
                                        26.5 0.1 .2 5 10 0.015 1.
  }
  2 {
    nDMaterial PressureIndependMultiYield 2 2 1.8 4.e4 2.e5 40 .1
  }
}

set gravY -$accGravity                 ;#calc. gravity
set gravX [expr -$gravY*$loadBias]

# define nodes
set numXnode  [expr $numXele+1]
set numYnode  [expr $numYele+1]

for {set i 1} {$i <= $numXnode} {incr i 1} {
  for {set j 1} {$j <= $numYnode} {incr j 1} {
    set xdim [expr ($i-1)*$xSize]
    set ydim [expr ($j-1)*$ySize] 
    set nodeNum [expr $i + ($j-1)*$numXnode] 
    node $nodeNum $xdim $ydim 
  }
}

# define elements
for {set i 1} {$i <= $numXele} {incr i 1} {
  for {set j 1} {$j <= $numYele} {incr j 1} {
    set eleNum [expr $i + ($j-1)*$numXele] 
    set n1  [expr $i + ($j-1)*$numXnode] 
    set n2  [expr $i + ($j-1)*$numXnode + 1] 
    set n4  [expr $i + $j*$numXnode + 1] 
    set n3  [expr $i + $j*$numXnode] 
                                    #      thick maTag  bulk  mDensity  perm1  perm2  gravity      press   
    element quadUP $eleNum $n1 $n2 $n4 $n3 1.0  $matOpt $bulk $fmass  $hperm  $vperm $gravX $gravY $press   
  }
}  

#set material to elastic for gravity loading
updateMaterialStage -material $matOpt -stage 0  

# fix the base, and free surface drainage
for {set i 1} {$i <= $numXnode} {incr i 1} {
  fix $i 1 1 0
  set surfnode [expr ($numYnode-1)*$numXnode + $i] 
  fix $surfnode 0 0 1
}

# tie all disp. DOFs at same level
for {set i 1} {$i < $numYnode} {incr i 1} {
  set nodeNum1 [expr $i*$numXnode + 1]
  for {set j 2} {$j <= $numXnode} {incr j 1} {
    set nodeNum2 [expr $i*$numXnode + $j]
    equalDOF $nodeNum1 $nodeNum2 1 2 
  }
}

#############################################################
# GRAVITY APPLICATION (elastic behavior)

# create the SOE, ConstraintHandler, Integrator, Algorithm and Numberer
numberer RCM
system ProfileSPD
test NormDispIncr 1.0e-8 25 2
algorithm Newton
constraints Penalty 1.e18 1.e18
integrator Newmark 1.5  1.  
analysis Transient 

analyze 3 5.e3

# switch material stage from elastic (gravity) to plastic
switch $matOpt {
  1 {
     updateMaterialStage -material $matOpt -stage 1
  }
  2 {
     updateMaterialStage -material $matOpt -stage 1
  }
}
analyze 5 5.e3
# rezero time
wipeAnalysis
setTime 0.0
#loadConst -time 0.0
#############################################################
# NOW APPLY LOADING SEQUENCE AND ANALYZE (plastic)

#   base input motion                     
pattern UniformExcitation  1  1    -accel "Sine 0. 10. $period -factor $accMul"

#input motion through data file
#pattern UniformExcitation  1  1  -accel "Series -factor $accMul -filePath $accNam -dt $accDt"

#recorder for nodal variables along the vertical center line.
set nodeList {}
for {set i 0} {$i < $numYnode} {incr i 1} {
  lappend nodeList [expr $numXnode/2 + $i*$numXnode]
}

#define recorders for disp., excess pore pressure, and acc.
#Note: disp and acc outputs are relative to the base
eval "recorder Node -file disp  -time -node $nodeList -dof 1 2 -dT $deltaT disp"
eval "recorder Node -file pwp  -time -node $nodeList -dof 3 -dT $deltaT vel"
eval "recorder Node -file acc  -time -node $nodeList -dof 1 2 -dT $deltaT accel"

#stress/strain output at Gauss point 1 of each element along center line
set name1 "stress";   set name2 "strain";  
for {set i 1} {$i < $numYnode} {incr i 1} {
  set ele [expr $numXele-$numXele/2+($i-1)*$numXele] 
  set name11 [join [list $name1 $i] {}]
  set name21 [join [list $name2 $i] {}] 
  recorder Element -ele $ele  -time -file $name11  -dT $deltaT material 1 stress
  recorder Element -ele $ele  -time -file $name21  -dT $deltaT material 1 strain
}

#analysis options
constraints Penalty 1.e18 1.e18
test NormDispIncr 1.e-5 25 0
numberer RCM
algorithm Newton
system ProfileSPD
#some mass proportional and initial-stiffness proportional damping
rayleigh $massProportionalDamping 0.0 $InitStiffnessProportionalDamping 0.0
integrator Newmark $gamma  [expr pow($gamma+0.5, 2)/4]  
analysis VariableTransient 

#analyze 
set startT [clock seconds]
analyze $numSteps $deltaT [expr $deltaT/100] $deltaT 15
set endT [clock seconds]
puts "Execution time: [expr $endT-$startT] seconds."

wipe  #flush ouput stream

MATLAB Plotting File

clear all;

a1=load('acc');
d1=load('disp');
p1=load('pwp');
s1=load('stress1');
e1=load('strain1');
s5=load('stress5');
e5=load('strain5');
s9=load('stress9');
e9=load('strain9');

fs=[0.5, 0.2, 4, 6];
accMul = 2;

%integration point 1 p-q
po=(s1(:,2)+s1(:,3)+s1(:,4))/3;
for i=1:size(s1,1)
    qo(i)=(s1(i,2)-s1(i,3))^2 + (s1(i,3)-s1(i,4))^2 +(s1(i,2)-s1(i,4))^2 + 6.0* s1(i,5)^2;
   qo(i)=sign(s1(i,5))*1/3.0*qo(i)^0.5;
end

figure(1); close 1; figure(1);
%integration point 1 stress-strain
subplot(2,1,1), plot(e1(:,4),s1(:,5),'r');
title ('shear stress \tau_x_y VS. shear strain \epsilon_x_y at 10 m depth');
xLabel('Shear strain \epsilon_x_y');
yLabel('Shear stress \tau_x_y (kPa)');
subplot(2,1,2), plot(-po,qo,'r');
title ('confinement p VS. deviatoric stress q at 10 m depth');
xLabel('confinement p (kPa)');
yLabel('q (kPa)');
set(gcf,'paperposition',fs);
saveas(gcf,'SS_PQ_10m','jpg');


%integration point 5 p-q
po=(s5(:,2)+s5(:,3)+s5(:,4))/3;
for i=1:size(s5,1)
    qo(i)=(s5(i,2)-s5(i,3))^2 + (s5(i,3)-s5(i,4))^2 +(s5(i,2)-s5(i,4))^2 + 6.0* s5(i,5)^2;
   qo(i)=sign(s5(i,5))*1/3.0*qo(i)^0.5;
end

figure(5); close 5; figure(5);
%integration point 5 stress-strain
subplot(2,1,1), plot(e5(:,4),s5(:,5),'r');
title ('shear stress \tau_x_y VS. shear strain \epsilon_x_y at 6 m depth');
xLabel('Shear strain \epsilon_x_y');
yLabel('Shear stress \tau_x_y (kPa)');
subplot(2,1,2), plot(-po,qo,'r');
title ('confinement p VS. deviatoric stress q at 6 m depth');
xLabel('confinement p (kPa)');
yLabel('q (kPa)');
set(gcf,'paperposition',fs);
saveas(gcf,'SS_PQ_6m','jpg');

%integration point 9 p-q
po=(s9(:,2)+s9(:,3)+s9(:,4))/3;
for i=1:size(s1,1)
    qo(i)=(s9(i,2)-s9(i,3))^2 + (s9(i,3)-s9(i,4))^2 +(s9(i,2)-s9(i,4))^2 + 6.0* s9(i,5)^2;
   qo(i)=sign(s9(i,5))*1/3.0*qo(i)^0.5;
end

figure(6); close 6; figure(6);
%integration point 9 stress-strain
subplot(2,1,1), plot(e9(:,4),s9(:,5),'r');
title ('shear stress \tau_x_y VS. shear strain \epsilon_x_y at 2 m depth');
xLabel('Shear strain \epsilon_x_y');
yLabel('Shear stress \tau_x_y (kPa)');
subplot(2,1,2), plot(-po,qo,'r');
title ('confinement p VS. deviatoric stress q at 2 m depth');
xLabel('confinement p (kPa)');
yLabel('q (kPa)');
set(gcf,'paperposition',fs);
saveas(gcf,'SS_PQ_2m','jpg');

figure(2); close 2; figure(2);
%node 3 displacement relative to node 1
subplot(2,1,1),a=plot(d1(:,1),d1(:,8),'r');
hold on
subplot(2,1,1),b=plot(d1(:,1),d1(:,14),'g');
hold on
subplot(2,1,1),c=plot(d1(:,1),d1(:,22),'b');
title ('Lateral displacement wrt base');
xLabel('Time (s)');
yLabel('Displacement (m)');
legend([a,b,c],'8m depth','4m depth', 'Surface',2)
set(gcf,'paperposition',fs);
saveas(gcf,'Disp','jpg');

s=accMul*sin(0:pi/50:20*pi);
s=[s';zeros(3000,1)];
s1=interp1(0:0.01:40,s,a1(:,1));

figure(3); close 3; figure(3);
%node acceleration
subplot(3,1,1),a=plot(a1(:,1),s1+a1(:,22),'r');
legend(a,'at surface',4);
title ('Lateral acceleration');
xLabel('Time (s)');
yLabel('Acceleration (m/s^2)');
subplot(3,1,2),a=plot(a1(:,1),s1+a1(:,14),'r');
legend(a,'4 m depth',4);
xLabel('Time (s)');
yLabel('Acceleration (m/s^2)');
subplot(3,1,3),a=plot(a1(:,1),s1+a1(:,8),'r');
legend(a,'8 m depth',4);
xLabel('Time (s)');
yLabel('Acceleration (m/s^2)');
set(gcf,'paperposition',fs);
saveas(gcf,'Acc','jpg');

figure(4); close 4; figure(4);
subplot(3,1,1),a=plot(p1(:,1),p1(:,11),'r');
legend(a,'1 m depth',4);
title ('Pore pressure');
xLabel('Time (s)');
yLabel('Pore pressure (kPa)');
subplot(3,1,2),a=plot(p1(:,1),p1(:,7),'r');
legend(a,'5 m depth',4);
xLabel('Time (s)');
yLabel('Pressure (kPa)');
subplot(3,1,3),a=plot(p1(:,1),p1(:,2),'r');
legend(a,'10 m depth',4);
xLabel('Time (s)');
yLabel('Pressure (kPa)');
set(gcf,'paperposition',fs);
saveas(gcf,'EPWP','jpg');