This is the final notebook for the regeneration of the simulation results for our SPE RSS paper, which is now submitted to the journal of natural gas (or sth like that). I have not included the gravity effect in the code, so only the results from cases 1 to 4 can be regenerated here.
In [ ]:
using JFVM, Roots, PyPlot, ProgressMeter, DataFrames, Dierckx, JLD
In [2]:
include("frac_flow_funcs.jl")
include("foam_flow_two_phase.jl")
Out[2]:
In [3]:
data1=readtable("TableA1.csv")
Out[3]:
In [84]:
case_no=1
fmmob=data1[:fmmob][case_no]
epdry=data1[:epdry][case_no]
fmdry=data1[:fmdry][case_no]
krw0=data1[:krw0][case_no]
krg0=data1[:krg0][case_no]
swc=data1[:swc][case_no]
sgc=data1[:sgc][case_no]
nw=data1[:nw][case_no]
ng=data1[:ng][case_no]
muw=data1[:muw][case_no]
mug=data1[:mug][case_no]
Out[84]:
In [85]:
# load Corey-type and foam relperms
(krw, krg, dkrwdsw, dkrgdsw, krgf, dkrgfdsw)=corey_rel_perms(
swc=swc, sgr=sgc, krg0=krg0, ng=ng, krw0=krw0,
nw=nw, fmmob=fmmob, epdry=epdry, fmdry=fmdry)
# Define the fractional flow functions
fw(sw)=(krg(sw)/mug)./(krg(sw)/mug+krw(sw)/muw)
fw_foam(sw)=(krgf(sw)/mug)./(krgf(sw)/mug+krw(sw)/muw)
# define the mobility function
labda_tot(sw)=krg(sw)/mug+krw(sw)/muw
labda_tot_foam(sw)=krgf(sw)/mug+krw(sw)/muw
Out[85]:
In [86]:
(xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj)=frac_flow_wag(muw=muw,
mug=mug, ut=1e-5, phi=0.2,
k=1e-12, swc=swc, sgr=sgc, krg0=krg0, ng=ng, krw0=krw0,
nw=nw, sw0=1.0, sw_inj=0.0, L=1, pv_inj=5)
Out[86]:
In [87]:
L=1.0
sw_int = Spline1D(xt_prf, sw_prf, k=1);
sw_int(0.0003)
Out[87]:
In [88]:
plot(t, p_inj)
xlabel("time [s]")
ylabel("Injection pressure [Pa]")
Out[88]:
In [90]:
plot(xt_prf, sw_prf, "o--")
ylabel("Water saturation")
xlabel("x/t [m/s]")
title("WAG", fontsize=18)
Out[90]:
In [91]:
(xt_shock, sw_shock, xt_prf, sw_prf, t, p_inj)=frac_flow_sag(
fmmob=fmmob, epdry=epdry, fmdry=fmdry,
muw=muw, mug=mug, ut=1e-5, phi=0.2,
k=1e-12, swc=swc, sgr=sgc, krg0=krg0, ng=ng, krw0=krw0,
nw=nw, sw0=1.0, sw_inj=0.0, L=1, pv_inj=5)
Out[91]:
In [93]:
plot(xt_prf, sw_prf, "o--")
ylabel("Water saturation")
xlabel("x/t [m/s]")
title("SAG", fontsize=18)
Out[93]:
In [94]:
plot(xt_prf, labda_tot_foam(sw_prf), "o--")
xlabel("x/t [m/s]")
ylabel("Total relative mobility [1/(Pa.s)]")
Out[94]:
In [95]:
plot(t, p_inj)
xlabel("time [s]")
ylabel("Injection pressure [Pa]")
Out[95]:
In [66]:
include("foam_flow_two_phase.jl")
Out[66]:
In [72]:
Nx = 80 # number of cells in x direction
Ny = 80 # number of cells in y direction
W = 1 # width
H = 1 # height
x=[linspace(0,0.5,70); linspace(0.51, 1.0, 10)]
y=collect(linspace(0,H,Ny))
m = createMesh1D(x)
m2=createMesh2D(Nx,Ny, W, H);
perm_ave=1e-12
perm_field=permfieldlogrnde(Nx, Ny, perm_ave, 0.1, 0.1, 0.1);
In [109]:
sw2=foam_stars(m2, perm_ave=perm_field, poros_ave=0.2, sw_init=sw_shock,
fmmob=fmmob, epdry=epdry, fmdry=fmdry,
muw=muw, mug=mug, swc=swc, sgr=sgc, krg0=krg0, ng=ng, krw0=krw0, nw=nw,
t_final=50, dt0=50.0/100);
In [111]:
#save("case1_shocksat_2d.jld", "sw", sw2)
visualizeCells(sw2)
colorbar()
title("water saturation, sw0=sw_shock")
Out[111]:
In [105]:
#save("case1_shocksat_2d.jld", "sw", sw2)
pcolor(labda_tot_foam(sw2.value[2:end-1,2:end-1]'), vmax=4000)
colorbar()
title("total relative mobility, sw0=sw_shock")
Out[105]:
In [ ]:
sw2_sat=foam_stars(m2, perm_ave=perm_field, poros_ave=0.2, sw_init=1.0,
fmmob=fmmob, epdry=epdry, fmdry=fmdry,
muw=muw, mug=mug, swc=swc, sgr=sgc, krg0=krg0, ng=ng, krw0=krw0, nw=nw,
t_final=50, dt0=50.0/500);
In [108]:
#save("case1_fullsat_2d.jld", "sw", sw2_sat)
visualizeCells(sw2_sat)
colorbar()
title("water saturation, sw0=1.0")
Out[108]:
In [107]:
# plot total relative mobility
pcolor(labda_tot_foam(sw2_sat.value[2:end-1,2:end-1]'))
colorbar()
title("total relative mobility, sw0=1.0")
Out[107]:
In [112]:
sw1=foam_stars(m, perm_ave=1e-12, poros_ave=0.2, sw_init=1.0,
fmmob=fmmob, epdry=epdry, fmdry=fmdry,
muw=muw, mug=mug, swc=swc, sgr=sgc, krg0=krg0, ng=ng, krw0=krw0,
nw=nw, t_final=50, dt0=50.0/100);
visualizeCells(sw1)
xlabel("x [m]")
ylabel("water saturation [-]")
Out[112]: