It is hard for us to type your scripts line by line because it is in figure format.
Here is another example we prepared for your reference:
model new
;---------------------------------------------------------------------
; One-dimensional earthquake excitation of uniform layer
;---------------------------------------------------------------------
fish automatic-create off
model large-strain off
model configure dynamic
model gravity 10
; Create zones
zone create brick size 1 1 20
;
fish define attach_gp
loop local ii (1,20)
local gp1 = gp.near(0,0,ii)
local gp2 = gp.near(1,0,ii)
local gp3 = gp.near(1,1,ii)
local gp4 = gp.near(0,1,ii)
local a1 = attach.create(gp2, gp1)
local a2 = attach.create(gp3, gp1)
local a3 = attach.create(gp4, gp1)
attach.snap(a1) = 0
attach.snap(a2) = 0
attach.snap(a3) = 0
endloop
end
[attach_gp]
;
zone cmodel assign plastic-hardening
zone property density 1000
zone property pressure-ref 100e3 o-c-r 1.2
zone property friction 30 dilation 6 cohesion 15e3
zone property stiffness-50-ref 8e6 stiffness-oed-ref 6e6 stiffness-ur-ref 24e6
zone property flag-small on stiffness-0-ref 120e6 strain-70 5e-4
;
zone initialize-stresses ratio 0.5
;
fish operator iniprin(z, modelname)
local pp = zone.pp(z)
if zone.model(z) = modelname then
local prin = zone.stress.prin(z)
zone.prop(z,‘stress-1-effective’) = prin->x + pp
zone.prop(z,‘stress-2-effective’) = prin->y + pp
zone.prop(z,‘stress-3-effective’) = prin->z + pp
endif
end
[iniprin(::zone.list,‘plastic-hardening’)]
; Boundary conditions
zone gridpoint fix velocity-y
zone gridpoint fix velocity-z
table ‘acc’ import ‘gilroy1.acc’
zone face apply acceleration-x -0.02 table ‘acc’ range position-z=0
; Histories
history interval 1
model history name=‘time’ dynamic time-total
zone history name=‘sxz1’ stress-xz zoneid = 1
zone history name=‘exz1’ strain-increment-xz zoneid = 1
zone history name=‘sxz10’ stress-xz zoneid = 10
zone history name=‘exz10’ strain-increment-xz zoneid = 10
zone history name=‘acc00’ acceleration-x position (0,0, 0)
zone history name=‘acc10’ acceleration-x position (0,0,10)
zone history name=‘acc20’ acceleration-x position (0,0,20)
;
model dynamic timestep fix 5e-4
model solve time-total 25
model save ‘EarthquakeExcitation-ph’
;model new
model large-strain on
; Single case of SeveralCyclicStrainLevels example
;
; Create zone
zone create brick size 1 1 1
; Assign model and properties
zone cmodel assign plastic-hardening
zone property density 1000 stiffness-50-reference=3.0e4 stiffness-ur-reference=8.0e4
zone property pressure-reference=100.0 exponent=0.55 poisson=0.25
zone property coefficient-normally-consolidation=0.40
zone property friction=40.0 dilation=10.0 cohesion=0.0
zone property stress-1-effective=-100.0 stress-2-effective=-100.0 …
stress-3-effective=-100.0
zone property flag-smallstrain=true stiffness-0-reference=30.0e4 …
strain-70=2e-4
; Boundary conditions
zone gridpoint fix velocity
zone face apply stress-xx=-100.0 range union position-x 0 position-x 1
zone face apply stress-yy=-100.0 range union position-y 0 position-y 1
zone initialize stress xx -100.0 yy -100.0 zz -100.0
;model solve
zone gridpoint initialize velocity-x [CycStrain/10.0] range position-z=1
; Dynamic settings
model dynamic timestep fix 1.0e-4
; Histories
zone history name=‘stress’ stress-xz position (0.5,0.5,0.5)
zone history name=‘strain’ strain-increment-xz position (0.5,0.5,0.5)
history interval 1
; Cyclic loading
model cycle 1000
zone gridpoint initialize velocity-x -1 multiply
model cycle 2000
zone gridpoint initialize velocity-x -1 multiply
model cycle 2000
; Export histories of result to table
history export ‘stress’ vs ‘strain’ table [‘PH-’+string(CaseNumber)]
; Clear model for next run
zone delete
history delete
This is for sig-3 model
;model new
model large-strain on
; Single case of SeveralCyclicStrainLevels example
;
; Create zone
zone create brick size 1 1 1
; Assign model and properties
zone cmodel assign elastic
zone property density 1000 shear 50000000 bulk 66666666.6666667
; Boundary conditions
zone gridpoint fix velocity
zone face apply stress-xx=-100.0 range union position-x 0 position-x 1
zone face apply stress-yy=-100.0 range union position-y 0 position-y 1
zone initialize stress xx -100.0 yy -100.0 zz -100.0
;model solve
zone gridpoint initialize velocity-x [CycStrain/10.0] range position-z=1
; Dynamic settings
model dynamic timestep fix 1.0e-4
zone dynamic damping hysteretic sigmoidal-3 1.0324 -0.6350 -1.4571 ...
reduction-minimum 0.02
; Histories
zone history name='stress' stress-xz position (0.5,0.5,0.5)
zone history name='strain' strain-increment-xz position (0.5,0.5,0.5)
history interval 1
; Cyclic loading
model cycle 1000
zone gridpoint initialize velocity-x -1 multiply
model cycle 2000
zone gridpoint initialize velocity-x -1 multiply
model cycle 2000
; Export histories of result to table
history export 'stress' vs 'strain' table ['sig3-'+string(CaseNumber)]
; Clear model for next run
zone delete
history delete
The above two models have same material parameters (Gmax and reference strain). In these models, Seed and Idrisss mean value curve is adopted to find reference strains.
This is the comparison of shear modulus reduction curves and damping curves with shear strain amplitude.
In this example, the PH model begins under isotropic stress, meaning the initial mobilized friction angle is zero and will increase over time. The model exhibits hardening with each successive cycle, causing unclosed loops. Additionally, since it uses the Mohr-Coulomb failure criterion, it shows asymmetric strengths in opposite directions. And, the PH model yields in shear failure during each cycle, leading to what appears to be a larger damping effect if we directly derive the damping from this plastic behavior, similar to the Mohr-Coulomb model. In sum, all these make the comparison to be less like apple to apple.
BTW, I am not sure if the method to estimate the damping based on the area ratio is still valid for this case of unclosed cycles.