How to iteratively model rock specimens with different joint orientations

Hi,
I want to create rock specimens with different joint orientations in 3DEC 9.0 to test the uniaxial compressive strength of rock containing various joint orientations. However, when I run my code, I get the error message: ‘Cannot cut blocks after cycling.’ How can I resolve this issue? Below is my code.
Thank you in advance!
JLUJYP

model new
block tolerance 1e-6
model random 10001
model mechanical timestep fix 0.005
model large-strain on
[global Size=1]
[global R=0.025Size]
[global High=0.1
Size]
block create tunnel radius=[R] length=0,[High] radius-ratio=20 dip=-90 dip-direction=0 blocks-radial=1 blocks-tangential=2 blocks-axial=1
block delete range position-x [-1R],[R] not
block delete range position-y [-1
R],[R] not
block join on
[global area = math.piR^2]
[global length = High]
;call ‘fish-support.fis’
fish define findList
gp_list_top = list
gp_list_bottom = list
loop foreach gp block.gp.list
if block.gp.isgroup(gp,‘top’)
gp_list_top(“end”) = gp
endif
if block.gp.isgroup(gp,‘bottom’)
gp_list_bottom(“end”) = gp
endif
endloop
end
fish define zzstress
findList
global gp_top = gp_list_top(1)
; compression positive,unit MPa
local top_reaction = -list.sum(block.gp.force.reaction.z(::gp_list_top))
local bottom_reaction = list.sum(block.gp.force.reaction.z(::gp_list_bottom))
zzstress = 0.5
(top_reaction + bottom_reaction)/area
end
; Function that determines if solving should stop
fish define halt
global gpHalt = block.gp.near(0,0,High)
DISP_threshold = High1.125e-4
DISP = math.abs(block.gp.disp.z(gpHalt ))
halt = DISP > DISP_threshold
end
model save ‘model_initial’
fish define solveAll
loop local beta (0,90,5) ; Check angles from 0-90 in 5 degree increments
command
block delete
block create tunnel radius=[R] length=0,[High] radius-ratio=20 dip=-90 dip-direction=0 blocks-radial=1 blocks-tangential=2 blocks-axial=1
block delete range position-x [-1
R],[R] not
block delete range position-y [-1R],[R] not
block join on
block cut joint-set dip @beta dip-direction 0 origin 0 0 [0.5
High]
block zone generate edgelength [0.5R]
block gridpoint group ‘top’ range position-z [High]
block gridpoint group ‘bottom’ range position-z 0
; Assign model and properties
block zone cmodel assign mohr-coulomb
block zone property bulk 1.e8 shea 7.e7 cohesion 2.e3 density 0.002
block zone property friction 40. dilation 0. tension 2400.
block contact jmodel mohr
block contact property stiffness-normal 7.6e7 stiffness-shear 7.6e7 cohesion 1.e3 tension 2000 friction 20
block contact material-table default property stiffness-normal 7.6e7 stiffness-shear 7.6e7 cohesion 1.e3 tension 2000 friction 20
; Assign boundary conditions
block gridpoint apply velocity-z [-1.25e-7
High] range group ‘top’
block gridpoint apply velocity-z 0 range group ‘bottom’
; Cycle till the target strain is reached
model solve fish-halt @halt
; Add results to table
table ‘result’ add (@beta,@zzstress)
end_command
end_loop
end
@solveAll
; Save the last state, and the accumulated table
model save ‘final’

As the error states, you cannot create further block cuts after you have initiated cycling?

In the above provided code, you create blocks and solve your initial model etc. (up to line 43). After your initial model, within your solveAll function you then delete all blocks and create new blocks for your “local beta” loop, with different joint orientations. This is all within a single model, without a model new command, hence I think the error.

To test sensitivity to joint orientations like this, the most efficient way would be to list your different sets of input parameters in an input file and use python to loop through that input file, passing the input parameters to 3DEC for iterative model runs. You’d then be able to restore all of your models, or saved outputs of different joint orientations to compare results.

Dear jmswdmn,

Thank you very much for your suggestions! I tried using ‘model new’ inside ‘solveAll’, but the software indicated that this is not allowed.

I am not sure how to use Python to loop through input files and run iterative modeling in 3DEC. Could you provide a simple example demonstration?

Thanks again for your help!

JLUJYP

simple example .py below, will load an input .csv file with variables x, y and z as a python dictionary. then for each row in the .csv file, pass these variables to a model.dat file to run.

caveat: not tested, but should get you 99% of the way there.

import itasca as it
it.command('python-reset-state false')
import csv
import io

var_dict = {}

with io.open("inputs.csv", 'r', encoding='utf-8-sig') as csvfile:
    rows = csv.DictReader(csvfile)
    col_to_var = {}
    for row_i, row in enumerate(rows):
        var_dict[row_i] = row

#for key, data in var_dict.items():
keys = list(var_dict.keys())
for key in keys[:]:
    print('working on %s' %key)
    x = float(data['x'])
    y = float(data['y'])
    z = float(data['z'])
    try:
        # Setup the 3DEC input parameters and call the model file to run
        it.command("""
        model new

[_x = {x}]
[_y = {y}]
[_z = {z}]

call "model.dat"
""".format(x = x, y = y, z = z))
    except RuntimeError:
        # If an error occurs in 3DEC, skip to next model
        print("Row %i failed to converge" %key)

        break

1 Like

Dear jmswdmn,
Thank you very much for your help! I will try your code.
Have a good day!
JLUJYP