Not Getting Expected Unloading Results PDMY Material

[armkam]

Hi,

I have modeled a simple medium sand soil column with PDMY material using four 1x1m quad elements which is fixed at based and constrained at 1st Dof (horizontally) at both sides.
A series of loading/unloading vertical load is applied on upper surface nodes. This is like compression test (Oedometer) and the loading phase of vertical stress-strain is expected but the unloading phase is so close to loading case which is so unexpected. The vertical load/unload is with 50KN increment increasing from 0 to 1000MN and then decreasing to 0 again.
The ideal curve and the obtained curve are as follows respectively:
https://i.imgsafe.org/d3b702703a.jpg

https://i.imgsafe.org/d45d09aa66.jpg

Also you can find the loading/unloading file if needed here:
http://www.fileconvoy.com/dfl.php?id=gf ... 02b213218b


And the script is copied here:
-------------------------------------------------------------------------------------------------------------------------------------------------
wipeAnalysis

model BasicBuilder -ndm 2 -ndf 2

node 1 0.0 0.0
node 2 1.0 0.0
node 3 0.0 1.0
node 4 1.0 1.0
node 5 0.0 2.0
node 6 1.0 2.0
node 7 0.0 3.0
node 8 1.0 3.0
node 9 0.0 4.0
node 10 1.0 4.0

fix 1 1 1
fix 2 1 1
fix 3 1 0
fix 4 1 0
fix 5 1 0
fix 6 1 0
fix 7 1 0
fix 8 1 0
fix 9 1 0
fix 10 1 0

set dt 5.0
set gamma 0.5

#---MATERIAL PROPERTIES
# soil mass density (Mg/m^3)
set rho [expr 18/9.81]
# soil shear velocity
set Vs 200.0
# soil elastic modulus
set Es 40000
# poisson's ratio of soil
set nu 0.3
# soil shear modulus
set Gmax [expr $Es/(2*(1+$nu))]
# soil bulk modulus
set bulk [expr $Es/(3*(1-2*$nu))]
# soil friction angle
set phi 35.0
# soil cohesion
set coh 0.0
# peak shear strain
set gammaPeak 0.1
# reference pressure
set refPress 80.0
# pressure dependency coefficient
set pressDependCoe 0.5
# phase transformation angle
set PTAng 27
# rate of shear-induced volume decrease (contraction) or pore pressure buildup
set contrac 0.07
# the rate of shear-induced volume increase (dilation)
set dilat1 0.4
set dilat2 2.0
# parameters controlling the mechanism of cyclic mobility (Kpa)
set liquefac1 10
set liquefac2 0.01
set liquefac3 1


nDMaterial PressureDependMultiYield 10 2 $rho $Gmax $bulk $phi $gammaPeak $refPress $pressDependCoe $PTAng $contrac $dilat1 $dilat2 $liquefac1 $liquefac2 $liquefac3
nDMaterial InitialStateAnalysisWrapper 1 10 2

set wgtX 0.0
set wgtY [expr -9.81*$rho]
element quad 1 1 2 4 3 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY
element quad 2 3 4 6 5 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY
element quad 3 5 6 8 7 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY
element quad 4 7 8 10 9 1.0 "PlaneStrain" 1 0.0 0.0 $wgtX $wgtY

recorder Node -file ESoilDispDir1.out -time -nodeRange 1 10 -dof 1 disp
recorder Node -file ESoilDispDir2.out -time -nodeRange 1 10 -dof 2 disp

recorder Node -file ESoilForceDir1.out -time -nodeRange 1 10 -dof 1 reaction
recorder Node -file ESoilForceDir2.out -time -nodeRange 1 10 -dof 2 reaction

recorder Element -file Estress1.out -time -dT $dt -eleRange 1 4 material 1 stress
recorder Element -file Estress2.out -time -dT $dt -eleRange 1 4 material 2 stress
recorder Element -file Estress3.out -time -dT $dt -eleRange 1 4 material 3 stress
recorder Element -file Estress4.out -time -dT $dt -eleRange 1 4 material 4 stress

recorder Element -file Estrain1.out -time -dT $dt -eleRange 1 4 material 1 strain
recorder Element -file Estrain2.out -time -dT $dt -eleRange 1 4 material 2 strain
recorder Element -file Estrain3.out -time -dT $dt -eleRange 1 4 material 3 strain
recorder Element -file Estrain4.out -time -dT $dt -eleRange 1 4 material 4 strain


constraints Transformation
test NormDispIncr 1.e-5 15 1
algorithm Newton
numberer RCM
system SparseGeneral
integrator LoadControl 1
analysis Static

InitialStateAnalysis on
updateMaterialStage -material 1 -stage 0

analyze 4

updateMaterialStage -material 1 -stage 1

analyze 2

InitialStateAnalysis off

setTime 0.0
wipeAnalysis

set Vloading Vloading.out
timeSeries Path 2 -dt $dt -filePath $Vloading -factor 0.5

pattern Plain 1 2 {
load 9 0.0 -1.0
load 10 0.0 -1.0
}

#---RAYLEIGH DAMPING PARAMETERS
set pi 3.141592654
# damping ratio
set damp 0.02
# lower frequency
set omega1 [expr 2*$pi*0.2]
# upper frequency
set omega2 [expr 2*$pi*20]
# damping coefficients
set a0 [expr 2*$damp*$omega1*$omega2/($omega1 + $omega2)]
set a1 [expr 2*$damp/($omega1 + $omega2)]
puts "damping coefficients: a_0 = $a0; a_1 = $a1"

recorder Node -file SoilDispDir1.out -time -nodeRange 1 10 -dof 1 disp
recorder Node -file SoilDispDir2.out -time -nodeRange 1 10 -dof 2 disp
recorder Node -file SoilSettlement.out -time -nodeRange 9 10 -dof 2 disp

recorder Node -file SoilForceDir1.out -time -nodeRange 1 10 -dof 1 reaction
recorder Node -file SoilForceDir2.out -time -nodeRange 1 10 -dof 2 reaction

recorder Element -file stress1.out -time -dT $dt -eleRange 1 4 material 1 stress
recorder Element -file stress2.out -time -dT $dt -eleRange 1 4 material 2 stress
recorder Element -file stress3.out -time -dT $dt -eleRange 1 4 material 3 stress
recorder Element -file stress4.out -time -dT $dt -eleRange 1 4 material 4 stress
recorder Element -file stressBase.out -time -dT $dt -eleRange 1 1 material 1 stress


recorder Element -file strain1.out -time -dT $dt -eleRange 1 4 material 1 strain
recorder Element -file strain2.out -time -dT $dt -eleRange 1 4 material 2 strain
recorder Element -file strain3.out -time -dT $dt -eleRange 1 4 material 3 strain
recorder Element -file strain4.out -time -dT $dt -eleRange 1 4 material 4 strain

constraints Penalty 1.e16 1.e16
test NormDispIncr 1.e-5 35 1
algorithm KrylovNewton
numberer RCM
system ProfileSPD
integrator Newmark $gamma [expr pow($gamma+0.5, 2.)/4.]
rayleigh $a0 $a1 0.0 0.0
analysis Transient
analyze 40010 $dt

[armkam]

I don't know whom I should ask this question and I have emailed code developers but got no answers so far.