Joint HM coupling: pore pressure diffusion stuck

3DEC
1

Hi there everyone!

So, we are running some tests for possible fault reactivation under given imposed pore pressure.

We started from the documentation example coupled-fast. When changing the model dimensions (from metric to kilometric scale), pressure does not propagate:

image
image1070×809 28.3 KB

We tried increasing the timestep and reducing the mesh size, but nothing changed.

Also, since we had to increase the mesh size, the apparent water bulk would greatly be reduced, hence getting out of the “incompressible fluid” necessity for fast-flow computations. We also ran the fully coupled script (without fast flow and with actual water bulk), but still no luck.

I provide the full script below.

Any advice welcome!

Cheers

Theophile

; fully coupled analysis with fast-flow
model new
model config fluid
model large-strain off

fish define set_props
    ;; Initial coupled-fast
    ;global xmin = -1.
    ;global xmax = 1.
    ;global ymin = xmin
    ;global ymax = xmax
    ;global zmin = xmin
    ;global zmax = xmax
    ;global z_jset = (zmin+zmax)/2.
    ;global edgelen = 0.25
    ;global ncycles = 2500
    
    ; Upscaled
    global xmin = -4000.
    global xmax = 4000.
    global ymin = -3000.
    global ymax = 3000.
    global zmin = -3000.
    global zmax = 0.
    global z_jset = (zmin+zmax)/2.
    global edgelen = 250.
    global ncycles = 10000

end
[set_props]

block create brick [xmin],[xmax] [ymin],[ymax] [zmin],[zmax]
block cut joint-set dip 0 dip-direction 90 origin 0 0 [z_jset]

block zone generate edgelength [edgelen]
block zone cmodel assign elastic
block zone prop density 2500 shear 2e9 bulk 5e9

block contact jmodel assign mohr
block contact property stiffness-normal 1e10 stiffness-shear 1e10 

flowplane vertex property aperture-initial 1e-4 ...
                          aperture-minimum 1e-5 aperture-maximum 1e-3

block fluid property bulk 2e9 density 1000 viscosity 1e-3

block gridpoint apply vel-x 0.0 range pos-x [xmin]
block gridpoint apply vel-x 0.0 range pos-x [xmax]
block gridpoint apply vel-y 0.0 range pos-y [ymin]
block gridpoint apply vel-y 0.0 range pos-y [ymax]
block gridpoint apply vel-z 0.0 range pos-z [zmin]
block gridpoint apply vel-z 0.0 range pos-z [zmax]

block insitu stress -1e7 -1e7 -1e7 0 0 0 nodis
block insitu p-p 2e6

; initial equilibrium - mechanical
model fluid active off
block fluid property bulk 0.0
model hist mech unbal-max
model cyc 100

; inital equilibrium - fluid
model fluid active on
model mech active off
block fluid property bulk 2e9
model cyc 100

block gridpoint ini disp 0 0 0
block contact reset disp 

flowknot apply p-p 3.0e6 range pos-x [xmin]
flowknot apply p-p 2.0e6 range pos-x [xmax]

block fluid fast-flow active on

model mech active on
block fluid property bulk [2e9/edgelen^3]
model mechanical substep 2000 
model mechanical slave on

history delete
model his fluid time-total
flowknot his p-p pos 0 0 [z_jset]

model fluid timestep max 1e-3
model solve fluid cycles [ncycles] or mechanical ratio 1e-6

;model save "coupled-fast"
2

The problem here is that your model is about 1000x bigger than the example, but the max fluid timestep is still set to 1e-3. If you change this to 1 (or 10, or 100) you will see some action. It will still take some time though - this mesh is quite a bit denser than the one in the example.

3

Hi @jhazzard , and sorry for the late response!

We did try increasing the timestep, but we did not increase it enough to get some effects.
So thanks for the hint!

Still, when fixing a too large timestep, high instabilities appear resulting in strange pore pressure distributions.
The figure below illustrates the pore pressure distribution through the joint + the pore pressure evolution at joint mid flow knot (dashed points), using the following parameters:

    global xmin = -1000.
    global xmax = 1000.
    global ymin = xmin
    global ymax = xmax 
    global zmin = xmin
    global zmax = xmax 
    global z_jset = (zmin+zmax)/2.
    global edgelen = 150.
    global ncycles = 60000

image
image1158×905 64.7 KB

We would have two additional questions:

  • according to the documentation, the maximum fluid flow timestep depends on so called permeability factors (see equation (19)). Is there any way to access those factors through 3DEC, Fish or Python?
  • the fast flow options requires an incompressible fluid. Following the example, fluid bulk is set through Kw/(mesh edge)^3. When having mesh edges greater than 1, we strongly deviate from the incompressible case. Is there anyway to bypass this? By, e.g., changing the reference unit system (kilometric)? After some testing, violating this condition with high mesh edges did not significantly change the results

Thanks a lot for your lights!

Regards

Theophile

PS : Meanwhile, a colleague of mine is experiencing a kind of adaptative timestep approach that sounds promising, I will forward her feedback!

4

The unfortunate truth is that with a very stiff fluid, and simple tetrahedral zones, the fluid-mechanical coupling is not very stable. Your best bet is one of the following:

  1. Don’t use the fast flow logic. Use regular coupled approach with a lower fluid modulus. Using 2e7 gave me the following for your model (no need to specify a maximum timestep in this case)
  2. Do an uncoupled approach as described [here]. (Joint Fluid Flow — Itasca Software 9.4 documentation). You can alternate back and forth between fluid-only and solid-only to get something sensible

image
image1141×622 46.3 KB