Data Files for Sliding Wedge” Verification Problem

make_system.dat

model new
model domain extent -2. 2. condition destroy
rblock create box -1. 1.
fracture create dip 40. dip-dir 130. position 0,0,1 size 5 
fracture create dip 60. dip-dir 220. position 0,-0.25,1 size 5  
rblock cut fractures 

; Boundary conditions 
rblock fix velocity spin range id 6 not

; Model Parameters 
rblock attribute density 2000.0 damp 0.7
contact cmat default model linear ...
                   property kn 1.e7 ks 1.e7 fric 1.e20
; Gravity
model gravity 0,0,-10.

; Monitoring
model history name 'aratio' mechanical ratio-average

; Settle under gravity loading to equilibrium state
model solve ratio-average 1e-5
model save 'ini.sav'
program return
Top

analytical_fric.py

import itasca as it 
import math
import numpy as np
from vec import vec3 

it.command("""
    model new
    model restore \'ini.sav\'
    python-reset-state off
      """)

f1=it.dfn.fracture.find(1)
f2=it.dfn.fracture.find(2)


# find normals to the fractures 
n1=f1.normal()
n2=f2.normal()

#Normal orientated to the interior of the bloc
if f1.normal_z()<0:
    n1=-n1
if f2.normal_z()<0:
    n2=-n2

#Apply Goodman and Shi theory

#Active force vector
WW=348.074          #find inertial mass from rigid block python object ???
A=vec3((0,0,-WW))

#Calculating the sliding direction
n12=n1.cross(n2)
s12=n12/n12.mag()*np.sign(n12.dot(A))

#Normal forces 
N1=-(A.cross(n2)).dot(n12)/(n12.mag())**2
N2=-(A.cross(n1)).dot(-n12)/(n12.mag())**2
#Shear force
T12=A.dot(s12)

# Critical friction angle (if zero cohesion)
fric_theory=math.atan(T12/(N1+N2))*180./math.pi
print("Analytical critical friction angle is "+str(fric_theory)+" degrees")
Top

pfc_fric.py

import itasca as it 
import math

fric_limit1 = 0.0
fric_limit2 = 90.

while True:

    it.command("""
    model new
    model restore \'ini.sav\'
    rblock attribute displacement 0
    rblock attribute velocity 0.
    rblock history name \'zdisp\' displacement-z id 6
    """)

    cycle_max = it.cycle() + 10000
  
    fric_angle = (fric_limit1 + fric_limit2)/2.
    it.command("""
                  [fric_angle = {}]
                  contact property fric [math.tan(fric_angle*math.pi/180.)])
                  model solve ratio-average 1e-6 cycles-total 10000
                  fish define isStable
                    isStable = false
                    if rblock.num == 4
                        if mech.cycle < 10000
                            isStable = true
                        endif
                    endif   
                  end
               """.format(fric_angle))
    
    if it.fish.call_function("isStable") is True: #stable
        fric_limit2 = fric_angle
        print("Stable for "+str(fric_angle) +" degrees")
        it.command("model save \'last_stable.sav\'")
    
    else : #unstable
        fric_limit1 = fric_angle
        print("Unstable for "+str(fric_angle) +" degrees")
        it.command("model save \'last_unstable.sav\'")
        
    #compute interval size
    min_interval = (fric_limit2 - fric_limit1)/2.
    if min_interval<0.001:
        break    
    
print("Last stable friction angle interval is ["
      + str("%.3f" % fric_limit1) +";"
      + str("%.3f" % fric_limit2)+"] (degrees)")
_error=abs(fric_limit2-fric_theory)/fric_theory*100.0
print("The error relative to the theoretical solution (" 
      + str("%.3f" % fric_theory)
      + ") is "+str("%.3f" % _error)+"%")
Top