; Model to test the effects of the contact between an impermeable and a permeable layers. 
; Model 5 = Model 3 which has an increased bulk modulus of the permeable layer compared to Model 2
; Model 5 has a sloping boundary between the two layers, 
; and a decreased cohesion and tensile strength of the low-perm. layer.
tit
FLAC Modelling Input File
;
config extra 4 gw
grid  100   100 
gen 0.00e0,0.00e0 0.00e0,1.00e4 3.0e4,1.00e4 3.0e4,0.0e0 i=1,101 j=1,101
mark i=1
mark j=1
mark i = 101 
mark j = 101

;
***************************************
model mohr
water bulk 2e9 dens 1000
prop poro=0.3 ; to be elaborated below
;; 1 mech step = 1 fluid step; modify if necessary!!!
set nmech 1 ngw 10 flow on
***************************************

tab 1 0, 0, 0, 10000, 13500, 10000, 16500, 5000, 30000, 5000, 30000, 0

gen tab 1


;
; Properties for each sub-regions: 
pro dens= 2720  bulk=4.09e9 shear=2.45e9 reg=20, 20
pro coh=3.00e7 tens=1.5e7 reg=20, 20
pro fric=30 dil= 5  reg=20, 20
pro perm=1.00e-14 porosity=0.2 reg=20, 20

; Properties for each sub-regions: 
pro dens= 2720  bulk=4.09e10 shear=2.45e9 reg=75, 75	; note higher bulk compared to perm. layer
pro coh=1.00e7 tens=1.5e6 reg=75, 75			; note the decreased coh and tens compared to perm. layer
pro fric=30 dil= 5  reg=75, 75
pro perm=1.00e-18 porosity=0.2 reg=75, 75

plo hold perm blo

;
apply syy = -2.e8 j= 101 
*****
ini sxx=-4.66800e8 var=0, 2.668e8 i=1, 101 j=1,101
ini syy=-4.66800e8 var=0, 2.668e8 i=1, 101 j=1,101
ini szz=-4.66800e8 var=0, 2.668e8 i=1, 101 j=1,101
ini sxy=0e6


fix x i=1
fix x i= 101 
fix y j=1

;
def hydro_pp
loop i (1,igp)
loop j (1,jgp)
if j = jgp Then
;sss2 = 3.67e7
sss2 = 7.35e7
else
aaa = abs(y(i, j) - y(i, jgp))
;sss2 = 1000 * 9.8 * aaa+3.67e7
sss2 = 1000 * 9.8 * aaa+7.35 e7 ;for 2kb
end_if
gpp(i, j) = sss2
end_loop
end_loop
end
hydro_pp
*

;
ini pp=7.35e7 j=101
fix pp j=101

ini sat 1
fix sat
;
set gravity=9.8
;
hist 1 unbal nstep 1

set small
set mech on fl off
step 200
ini   xd 0    yd 0
ini xvel 0  yvel 0
;plot hold hist 1 
set mech on fl on
step 200  ; this can be increased if you think the unbalanced forces haven't bottomed out
ini   xd 0    yd 0 
ini xvel 0  yvel 0

set large
plot hold hist 1


;initial velocity command should be here ie ini xvel 0.004 etc;


;
set nmech 1 ngw 1 flow on
ini xvel  0.002 var=-0.004 ,0 i=1, 101  j=1, 101 
step 30000
save compl5_0001.sav
step 90000
save compl5_0002.sav

set nmech 1 ngw 10 flow on


ini xvel  0.0002 var=-0.0004 ,0 i=1, 101  j=1, 101 

;set mech off
apply dis= 1e-10 i= 1 101 j= 1


def files
    if nn<10 then
     savefile = 'compl5_00' + string(nn) + '.sav'
    else
    if nn>9 then
     if nn<100 then
     savefile = 'compl5_0' + string(nn) + '.sav'
    else
    if nn>99 then
     savefile = 'compl5_' + string(nn) + '.sav'
    end_if
    end_if
    end_if
    end_if
    COMMAND
        save @savefile
    END_COMMAND
end

define step_loop
loop nn (1,60)
   ;loop n (1,1000)
	command
 		step 1500
	end_command
      ;vsi_perm
   ;end_loop
files
end_loop
end
step_loop




return