PressureDependMultiYield02-Example 3
Jump to navigation
Jump to search
- Command_Manual
- Tcl Commands
- Modeling_Commands
- model
- uniaxialMaterial
- ndMaterial
- frictionModel
- section
- geometricTransf
- element
- node
- sp commands
- mp commands
- timeSeries
- pattern
- mass
- block commands
- region
- rayleigh
- Analysis Commands
- Output Commands
- Misc Commands
- DataBase Commands
Input File
# single 20_8_BrickUP element with pressure dependent material.
# subjected to 1D sinusoidal base shaking
# Written by Jinchi Lu (May 2004)
set matOpt 1 ;# 1 = pressure depend;
;# 2 = pressure independ;
wipe
set friction 31. ;#friction angle
set phaseTransform 26. ;#phase transformation angle
set G1 9.e4 ;
set B1 22.e4 ;
set gamma 0.600 ;# Newmark integration parameter
set dT 0.01 ;# time step for analysis, does not have to be the same as accDt.
set numSteps 2500 ;# number of time steps
set rhoS 1.80 ;# saturated mass density
set rhoF 1.00 ;# fluid mass density
set Bfluid 2.2e6 ;# fluid shear modulus
set fluid1 1 ;# fluid material tag
set solid1 10 ;# solid material tag
set perm 1.e-5 ;#permeability (m/s)
set accGravity 9.81 ;#acceleration of gravity
set perm [expr $perm/$accGravity/$rhoF] ;# actual value used in computation
set accMul 2 ;# acceleration multiplier
set pi 3.1415926535 ;
set inclination 0;
set massProportionalDamping 0.0 ;
set InitStiffnessProportionalDamping 0.003;
set gravityX [expr $accGravity*sin($inclination/180.0*$pi)] ;# gravity acceleration in X direction
set gravityY 0.0 ;# gravity acceleration in Y direction
set gravityZ [expr -$accGravity*cos($inclination/180.0*$pi)] ;# gravity acceleration in Z direction
set ndm 3 ;# space dimension
model BasicBuilder -ndm 3 -ndf 4
node 1 0.00000 0.0000 0.00000
node 2 1.00000 0.0000 0.00000
node 3 1.00000 1.0000 0.00000
node 4 0.00000 1.0000 0.00000
node 5 0.00000 0.0000 1.00000
node 6 1.00000 0.0000 1.00000
node 7 1.00000 1.0000 1.00000
node 8 0.00000 1.0000 1.00000
fix 1 1 1 1 0
fix 2 1 1 1 0
fix 3 1 1 1 0
fix 4 1 1 1 0
fix 5 0 1 0 1
fix 6 0 1 0 1
fix 7 0 1 0 1
fix 8 0 1 0 1
model BasicBuilder -ndm 3 -ndf 3
node 9 0.50000 0.0000 0.00000
node 10 1.00000 0.5000 0.00000
node 11 0.50000 1.0000 0.00000
node 12 0.00000 0.5000 0.00000
node 13 0.50000 0.0000 1.00000
node 14 1.00000 0.5000 1.00000
node 15 0.50000 1.0000 1.00000
node 16 0.00000 0.5000 1.00000
node 17 0.00000 0.0000 0.50000
node 18 1.00000 0.0000 0.50000
node 19 1.00000 1.0000 0.50000
node 20 0.00000 1.0000 0.50000
fix 9 1 1 1
fix 10 1 1 1
fix 11 1 1 1
fix 12 1 1 1
fix 13 0 1 0
fix 14 0 1 0
fix 15 0 1 0
fix 16 0 1 0
fix 17 0 1 0
fix 18 0 1 0
fix 19 0 1 0
fix 20 0 1 0
# equalDOF
# tied nodes around
equalDOF 5 6 1 3
equalDOF 5 7 1 3
equalDOF 5 8 1 3
equalDOF 5 13 1 3
equalDOF 5 14 1 3
equalDOF 5 15 1 3
equalDOF 5 16 1 3
equalDOF 17 18 1 3
equalDOF 17 19 1 3
equalDOF 17 20 1 3
# define material and properties
switch $matOpt {
1 {
nDMaterial PressureDependMultiYield02 1 3 $rhoS $G1 $B1 $friction .1 80 0.5\
$phaseTransform 0.067 0.23 0.06 0.27
}
2 {
nDMaterial PressureIndependMultiYield 2 3 1.8 4.e4 2.e5 40 .1
}
}
element 20_8_BrickUP 1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 $matOpt $Bfluid $rhoF $perm $perm $perm $gravityX $gravityY $gravityZ
#recorder for nodal variables along the vertical center line.
set SnodeList {}
for {set i 0} {$i < 20} {incr i 1} {
lappend SnodeList [expr $i+1]
}
set FnodeList {}
for {set i 0} {$i < 8} {incr i 1} {
lappend FnodeList [expr $i+1]
}
# GRAVITY APPLICATION (elastic behavior)
numberer Plain
system ProfileSPD
test NormDispIncr 1.0e-8 20 1
algorithm KrylovNewton
constraints Penalty 1.e18 1.e18; #
set nw 1.5
set nw2 [expr pow($nw+0.5, 2)/4]
integrator Newmark $nw $nw2
analysis Transient
analyze 10 5.e3
# switch the material to plastic
updateMaterialStage -material $matOpt -stage 1
analyze 10 1.e1
setTime 0.0 ;# reset time, otherwise reference time is not zero for time history analysis
wipeAnalysis
eval "recorder Node -file disp -time -node $SnodeList -dof 1 2 3 -dT $dT disp"
eval "recorder Node -file pwp -time -node $FnodeList -dof 4 -dT $dT vel"
eval "recorder Node -file acc -time -node $SnodeList -dof 1 2 3 -dT $dT accel"
recorder Element -ele 1 -time -file stress1 -dT $dT material 1 stress
recorder Element -ele 1 -time -file strain1 -dT $dT material 1 strain
recorder Element -ele 1 -time -file stress5 -dT $dT material 5 stress
recorder Element -ele 1 -time -file strain5 -dT $dT material 5 strain
recorder Element -ele 1 -time -file stress17 -dT $dT material 17 stress
recorder Element -ele 1 -time -file strain17 -dT $dT material 17 strain
############# create dynamic time history analysis ##################
pattern UniformExcitation 1 1 -accel "Sine 0 10 1 -factor $accMul"
rayleigh $massProportionalDamping 0.0 $InitStiffnessProportionalDamping 0.0
integrator Newmark $gamma [expr pow($gamma+0.5, 2)/4]
constraints Penalty 1.e18 1.e18 ;# can't combine with test NormUnbalance
test NormDispIncr 1.0e-3 25 0 ;# can't combine with constraints Lagrange
#algorithm Newton ;# tengent is updated at each iteration
algorithm KrylovNewton ;# step not each iteration
system ProfileSPD ;# Use sparse solver. Next numberer is better to be Plain.
numberer Plain ;# method to map between between equation numbers of DOFs
analysis VariableTransient ;# splitting time step requires VariableTransient
############# perform the Analysis and record time used #############
set startT [clock seconds]
analyze $numSteps $dT [expr $dT/64] $dT 15
set endT [clock seconds]
puts "Execution time: [expr $endT-$startT] seconds."
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('stress17');
e9=load('strain17');
fs=[0.5, 0.2, 4, 6];
fs2=[0.5, 0.2, 4, 3];
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 +s1(i,6)^2+s1(i,7)^2) ;
qo(i)=sign(s1(i,7))*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(:,7),s1(:,7),'r');
title ('shear stress \tau_x_z VS. shear strain \epsilon_x_z at integration point 1');
xLabel('Shear strain \epsilon_x_z');
yLabel('Shear stress \tau_x_z (kPa)');
subplot(2,1,2), plot(-po,qo,'r');
title ('confinement p VS. deviatoric stress q at integration point 1');
xLabel('confinement p (kPa)');
yLabel('q (kPa)');
set(gcf,'paperposition',fs);
saveas(gcf,'SS_PQ_p1','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 + s5(i,6)^2 + s5(i,7)^2);
qo(i)=sign(s5(i,7))*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(:,7),s5(:,7),'r');
title ('shear stress \tau_x_z VS. shear strain \epsilon_x_z at integration point 5');
xLabel('Shear strain \epsilon_x_z');
yLabel('Shear stress \tau_x_z (kPa)');
subplot(2,1,2), plot(-po,qo,'r');
title ('confinement p VS. deviatoric stress q at integration point 5');
xLabel('confinement p (kPa)');
yLabel('q (kPa)');
set(gcf,'paperposition',fs);
saveas(gcf,'SS_PQ_p5','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 + s9(i,6)^2 + s9(i,7)^2);
qo(i)=sign(s9(i,7))*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(:,7),s9(:,7),'r');
title ('shear stress \tau_x_z VS. shear strain \epsilon_x_z at integration point 17');
xLabel('Shear strain \epsilon_x_z');
yLabel('Shear stress \tau_x_z (kPa)');
subplot(2,1,2), plot(-po,qo,'r');
title ('confinement p VS. deviatoric stress q at integration point 17');
xLabel('confinement p (kPa)');
yLabel('q (kPa)');
set(gcf,'paperposition',fs);
saveas(gcf,'SS_PQ_p17','jpg');
figure(2); close 2; figure(2);
%node 3 displacement relative to node 1
plot(d1(:,1),d1(:,14));
title ('Lateral displacement at element top');
xLabel('Time (s)');
yLabel('Displacement (m)');
set(gcf,'paperposition',fs2);
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
a = plot(a1(:,1),s1+a1(:,14),'r');
title ('Lateral acceleration at element top');
xLabel('Time (s)');
yLabel('Acceleration (m/s^2)');
set(gcf,'paperposition',fs2);
saveas(gcf,'Acc','jpg');
figure(4); close 4; figure(4);
a=plot(p1(:,1),p1(:,2));
title ('Pore pressure at base');
xLabel('Time (s)');
yLabel('Pore pressure (kPa)');
set(gcf,'paperposition',fs2);
saveas(gcf,'EPWP','jpg');
Displacement Output File
Stress-Strain Output File (integration point 1)
Stress-Strain Output File (integration point 5)
Stress-Strain Output File (integration point 17)
Excess Pore Pressure Output File
Acceleration Output File
Return to:
- NDMaterial Command
- UC San Diego soil models (Linear/Nonlinear, dry/drained/undrained soil response under general 2D/3D static/cyclic loading conditions (please visit UCSD for examples)
- UC San Diego Saturated Undrained soil
- Element Command
- UC San Diego u-p element (saturated soil)
- Related References