(Apparent?) tensile rupture in joint with a hydrostatic model

3DEC
1

Hi everyone!

Doing some tests on hydromechanics in 3DEC, I came up ith this situation I cannot explain: when considering pore pressure in joints (fluid not configured), I get tensile rupture even when specifying extreme tensile strength (10^10 Pa versus 10^7 Pa pore pressure).

If I consider a model with fluid configured and flow off, the tensile ruptures disappear.

Documentation mentions that when model is not configured for fluid, pore pressure are only used to estimate elastoplastic state:

The user can simply specify water pressure everywhere in the model and these pressures are used to calculate effective stresses. Effective stresses are then used in the calculations to determine if there is failure of solid material or slip on joints.

However, I would not expect rupture to appear when pore pressure is way below normal stress and tensile strength.

Which brings my question: is the tensile rupture state apparent or actual (hence influencing mechanical equilibrium resolution)?

I provide a detailed exemple in PS of this post.

Thanks a lot!

Cheers

Theophile

PS :
I considered a really simple scenario:

  • 1x1x1 km^3 cube
  • single joint
  • elastic matrix and elastoplastic joint (Mohr-Coulomb)
  • gravity
  • oedometric boundary conditions (i.e., blocked normal velocity on all faces except top one which is free)
  • 3-step balancing (with displacements and state resetting between each):
    • elastic joint + no pore pressure
    • elastoplastic joint + no pore pressure
    • elastoplastic joint + hydrostatic gradient

image
image1333×742 292 KB

Normal stress seems to be affected (poorly, though).

image
image1329×744 209 KB

2

My guess is that the joint failed in shear first. This would cause the tensile strength to drop to zero. If this is not the answer, please post your simple data file.

3

Thank you, @jhazzard !

You guessed right: adding residual tensile strength preserves from having tensile rupture (while having only tensile strength alone did not).

We still get some differences in contours but these are acceptable.

image
image809×795 28.1 KB

image
image835×838 29 KB

I still do not understand why no tensile rupture appears when considering config fluid/flow off, even with no residual tensile strength?
Aren’t we supposed to be in the same conditions as no config fluid?

Just in case, I provide below the scripts for the study

Thanks!

Theophile

PS :

Script WITHOUT configure fluid 
; Testing pore pressure effects on tensile rupture: no config fluid
; Model
; - elastic matrix
; - elastoplastic joint
; - gravity + oedometric BC
; - hydrostatic profile
; - 3-phase equilibrium
;  o elastic, no pore pressure
;  o plastic, no pore pressure
;  o plastic, pore pressure
; Results
; --> tensile rupture appears except if one gives a residual tensile strength (suggestion by J. Hazzard)

project new
model large-strain off
block mechanical mass-scale on
block mechanical damping global

; GEOMETRY
; ;;;;;;;;

; Parameters
def geom_parameters
    global xmin=0
    global ymin=0
    global zmin=-1000
    global xmax=1000
    global ymax=1000
    global zmax=0
    global fault = map('origin', vector(500,500,-500), 'dip', 70, 'dip-dir', 270)
end
[geom_parameters]

; Ranges
; - main fault
model range create 'right' &
    plane origin [fault('origin')->x], [fault('origin')->y], [fault('origin')->z] &
    dip [fault('dip')] dip-direction [fault('dip-dir')] below
model range create 'left' named-range 'right' not

; Geometry+mesh
block create brick [xmin],[xmax] [ymin],[ymax] [zmin],[zmax]
block cut joint-set &
    origin [fault('origin')->x], [fault('origin')->y], [fault('origin')->z] &
    dip [fault('dip')] dip-direction [fault('dip-dir')]
block zone gen edge 50

; Groups
block group 'side=right' range named-range 'right'
block group 'side=left' range named-range 'left'
block contact group 'contact=fault' range group-intersection 'side=right' 'side=left'

; Boundaries
block face group 'west' range pos-x [xmin]
block face group 'east' range pos-x [xmax]
block face group 'south' range pos-y [ymin]
block face group 'north' range pos-y [ymax]
block face group 'bottom' range pos-z [zmin]
block gridpoint group 'east' range pos-x [xmax]


; RUN: PHASE 1/3 ELASTIC, NO FLUID
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
def hm_parameters
    ; Gravity
    global grav_z = 9.81
    
    ; Zones
    global young_ = 56.5e9
    global nu_ = 0.28
    global rho_ = 2680
    ; Contacts
    global kn_ = 10.e9
    global ks_ = 5.e9
    global cohe_ = 0.
    global phi_ = 15.
    global psi_ = 5.
    global usc_ = 1.
    ; NB: turn to non-zero fictitious value to avoid tensile rupture
    global rt_tres_ = 0
    
    ; Fluid
    global density_water = 1000
    global bulk_water = 0 ;3.8e5
    global viscosity_water = 1e-3
    global pp_top = 0 ; Pore pressure at model top
    global pp_grad = -math.abs(density_water*grav_z) ; increase with depth (z<0)
    global pp_z0 = pp_top - zmax*pp_grad ; pore pressure at z=0
end
[hm_parameters]

; - generic
model gravity [grav_z]
; - zones
block zone cmodel assign elastic
block zone property young [young_] poisson [nu_] dens [rho_]
; - contacts
block contact jmodel assign elastic
block contact property &
    stiffness-normal [kn_] &
    stiffness-shear [ks_]

; BC = oedometric
def sigH_on_sigV(nu)
    ; SigmaH / SigmaV ratio considering an overburden in an elastic, homogeneous
    ; and isotropic material
    ; nu: float, Poisson ratio
    sigH_on_sigV = nu/(1-nu)
end
block insitu topography ratio-x [sigH_on_sigV(nu_)] ratio-y [sigH_on_sigV(nu_)] ratio-z 1
block face apply velocity-normal 0.0 range group 'west'
block face apply velocity-normal 0.0 range group 'east'
block face apply velocity-normal 0.0 range group 'south'
block face apply velocity-normal 0.0 range group 'north'
block face apply velocity-normal 0.0 range group 'bottom'

model solve ratio 1.e-6
model save 'no_config_fluid_1-on-3.3dsav'

; RUN: PHASE 2/3 PLASTIC, NO FLUID
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
block gridpoint initialize displacement 0 0 0
block contact reset displacement
block contact jmodel mohr range group "contact=fault"
block contact property &
    stiffness-normal [kn_] &
    stiffness-shear [ks_] &
    cohesion [cohe_] &
    friction [phi_] &
    dilation [psi_] &
    dilation-zero = [usc_] &
    tension-residual = [rt_tres_] &
    range group "contact=fault"
model solve ratio 1.e-6
model save 'no_config_fluid_2-on-3.3dsav'

; RUN: PHASE 3/3 PLASTIC, FLUID
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
block gridpoint initialize displacement 0 0 0
block contact reset displacement
block contact reset state
; Hydrostatic
block insitu &
    pore-pressure [pp_z0] &
    gradient (0,0,[pp_grad])
model solve ratio 1.e-6
model save 'no_config_fluid_3-on-3.3dsav'
Script WITH configure fluid/flow off 
; Testing pore pressure effects on tensile rupture: with config fluid
; Model
; - elastic matrix
; - elastoplastic joint
; - gravity + oedometric BC
; - hydrostatic profile
; - 3-phase equilibrium
;  o elastic, no pore pressure
;  o plastic, no pore pressure
;  o plastic, pore pressure
; Results
; --> no tensile rupture

project new
model large-strain off
model configure fluid
block mechanical mass-scale on
block mechanical damping global

; GEOMETRY
; ;;;;;;;;

; Parameters
def geom_parameters
    global xmin=0
    global ymin=0
    global zmin=-1000
    global xmax=1000
    global ymax=1000
    global zmax=0
    global fault = map('origin', vector(500,500,-500), 'dip', 70, 'dip-dir', 270)
end
[geom_parameters]

; Ranges
; - main fault
model range create 'right' &
    plane origin [fault('origin')->x], [fault('origin')->y], [fault('origin')->z] &
    dip [fault('dip')] dip-direction [fault('dip-dir')] below
model range create 'left' named-range 'right' not

; Geometry+mesh
block create brick [xmin],[xmax] [ymin],[ymax] [zmin],[zmax]
block cut joint-set &
    origin [fault('origin')->x], [fault('origin')->y], [fault('origin')->z] &
    dip [fault('dip')] dip-direction [fault('dip-dir')]
block zone gen edge 50

; Groups
block group 'side=right' range named-range 'right'
block group 'side=left' range named-range 'left'
block contact group 'contact=fault' range group-intersection 'side=right' 'side=left'

; Boundaries
block face group 'west' range pos-x [xmin]
block face group 'east' range pos-x [xmax]
block face group 'south' range pos-y [ymin]
block face group 'north' range pos-y [ymax]
block face group 'bottom' range pos-z [zmin]
block gridpoint group 'east' range pos-x [xmax]


; RUN: PHASE 1/3 ELASTIC, NO FLUID
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
def hm_parameters
    ; Gravity
    global grav_z = 9.81
    
    ; Zones
    global young_ = 56.5e9
    global nu_ = 0.28
    global rho_ = 2680
    ; Contacts
    global kn_ = 10.e9
    global ks_ = 5.e9
    global cohe_ = 0.
    global phi_ = 15.
    global psi_ = 5.
    global usc_ = 1.
    
    ; Fluid
    global density_water = 1000
    global bulk_water = 0 ;3.8e5
    global viscosity_water = 1e-3
    global pp_top = 0 ; Pore pressure at model top
    global pp_grad = -math.abs(density_water*grav_z) ; increase with depth (z<0)
    global pp_z0 = pp_top - zmax*pp_grad ; pore pressure at z=0
end
[hm_parameters]

; - generic
model gravity [grav_z]
block fluid property bulk [bulk_water]
block fluid property density [density_water]
block fluid property viscosity [viscosity_water]
; - zones
block zone cmodel assign elastic
block zone property young [young_] poisson [nu_] dens [rho_]
; - contacts
block contact jmodel assign elastic
block contact property &
    stiffness-normal [kn_] &
    stiffness-shear [ks_]
flowplane vertex property &
    aperture-initial 1e-5 &
    aperture-minimum 1e-5 &
    aperture-maximum 1e-5
    
model fluid off mech on
model cycle 0

; BC = oedometric
def sigH_on_sigV(nu)
    ; SigmaH / SigmaV ratio considering an overburden in an elastic, homogeneous
    ; and isotropic material
    ; nu: float, Poisson ratio
    sigH_on_sigV = nu/(1-nu)
end
block insitu topography ratio-x [sigH_on_sigV(nu_)] ratio-y [sigH_on_sigV(nu_)] ratio-z 1
block face apply velocity-normal 0.0 range group 'west'
block face apply velocity-normal 0.0 range group 'east'
block face apply velocity-normal 0.0 range group 'south'
block face apply velocity-normal 0.0 range group 'north'
block face apply velocity-normal 0.0 range group 'bottom'

model solve ratio 1.e-6
model save 'config_fluid_flow_off_1-on-3.3dsav'

; RUN: PHASE 2/3 PLASTIC, NO FLUID
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
block gridpoint initialize displacement 0 0 0
block contact reset displacement
block contact jmodel mohr range group "contact=fault"
block contact property &
    stiffness-normal [kn_] &
    stiffness-shear [ks_] &
    cohesion [cohe_] &
    friction [phi_] &
    dilation [psi_] &
    dilation-zero = [usc_] &
    range group "contact=fault"
model solve ratio 1.e-6
model save 'config_fluid_flow_off_2-on-3.3dsav'

; RUN: PHASE 3/3 PLASTIC, FLUID
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
block gridpoint initialize displacement 0 0 0
block contact reset displacement
block contact reset state
; Hydrostatic
block insitu &
    pore-pressure [pp_z0] &
    gradient (0,0,[pp_grad])
model solve ratio 1.e-6
model save 'config_fluid_flow_off_3-on-3.3dsav'
4

OK - it seems I was not correct. Your tensile strength is 0, so the fact that the joint fails in shear is irrelevant. The real problem is the global damping. This should never be used - it is only still in the code for people who want to try to reproduce model results from the old days. In the model with no fluid and global damping, there is a bit of instability when the pore pressure changes, which causes a bunch of tensile failure. Just comment out line 17, and the results are much more similar.

5

It seems to me you were correct!
Turning residual tensile strength to non-zero does prevent from having tension rupture, even with a zero tensile strength.
It makes sense because, as you presumed, slip rupture happens first and 3DEC directly switches to residual strength.

Turning global damping off indeed brings very similar results to the model with config flow.
However, when turning residual tensile strength to 0 for the no config fluid model, there is still tension rupture.

image
image767×850 37.1 KB

To sum up - and taking the model with config fluid, no global damping and zero residual tensile strength as reference (this model has no tensile rupture) - we get for the no config fluid model:

  • with global damping and with zero residual tensile strength: tensile rupture; more or less comparable shear/normal results to reference
  • with global damping and with non-zero residual tensile strength: no tensile rupture; more or less comparable shear/normal results to reference
  • without global damping and with zero residual tensile strength: tensile rupture; very comparable shear/normal results to reference
  • without global damping and with non-zero residual tensile strength: no tensile rupture; very comparable shear/normal results to reference

The 4th model is the best candidate.

Theophile

6

Probably the timestep/mass density scaling is a bit different when you turn on fluid, so the abrupt application of pore pressure causes tensile failure in one case and not the other. I was running this in version 9, where the default damping is “combined” and the results were more similar.

7

Thanks a lot for your clarification, I’ll dig in this direction!

Cheers

Theophile