Verification of dynamic performance of PH-model model

hello. We are checking the dynamic performance of the PH-model.

I used it with flag-smallstrain on.
Also, the comparison model used cmodel elastic & sigmoidal-3.

Here, I want to check whether the result of this comparison is correct or incorrect.

Code and pictures attached.
thank you

image
image


In the graph, the Y-axis is stress and the X-axis is strain. Multiply X by 100 to get a %.

Also, flac3d ver. I used the 9.00 demo version.

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’

And some results:

Acc at surface

Stress-strain at base

Stress-strain in the middle

The two models (PH-small and sig3 model) are tested using the example code "One-Zone Sample Exercised at Several Cyclic Strain Levels’ ( One-Zone Sample Exercised at Several Cyclic Strain Levels — Itasca Software 9.0 documentation (itascacg.com))

The following code is for PH-small 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 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.

The PH cycles: