In this notebook, I'll read all the sheets in the excel files provide by Hamid, and save them in csv files. Then I'll load and analyze them.
For a weird reason, all the units in this notebook are SI, except for energy and exergy that are reported in kJ instead of J.
Sometimes, I forget to say that the units are kJ/s and not kJ. I try to fix, otherwise figure it out for yourselves.
BE CAREFUL: always run this notebook by choosing run all from the cell menu.
In [16]:
using JLD
# efficiency factors
η_pc = 0.34 # pulverized coal power plant efficiency
η_igcc = 0.40 # integrated gasification combined cycle
η_ngcc = 0.45 # natural gas combined cycle
# fossil fuels chemical exergy
MW_CO2 = 0.044 # kg/mol
# 1- coal; formula: CH
ex_coal = 400 + 100 # kJ/mol
MW_coal = 12 + 1 # g/mol
em_coal = MW_CO2/ex_coal # kg CO2/kJ chemical exergy
# 2- natural gas; formula: CH4
ex_gas = 400 + 4*100 # kJ/mol
MW_gas = 12 + 4 # g/mol
em_gas = MW_CO2/ex_gas # kg CO2/kJ chemical exergy
e_pc = em_coal/η_pc # kg CO2/kJ electricity
e_igcc = em_coal/η_igcc # kg CO2/kJ electricity
e_ngcc = em_gas/η_ngcc # kg CO2/kJ electricity
eta_pump = 0.8; # pump mechanical efficiency
eta_driver = 0.8 # electrical driver efficiency
ex_ccs = 5000 # [kJ/kg CO2]
Out[16]:
In [17]:
# @save "input_param.jld"
The above numbers will be used later in the exergy analysis below.
In [18]:
using DataFrames, PyPlot, CoolProp, Dierckx, CSV
include("SetPyPlot.jl")
SetPyPlot.setrcparam()
PyPlot.svg(false)
p0 = 101325 # [Pa]
T0 = 20+273.15 # [K]
Out[18]:
In [19]:
# some functions
function fric_factor(Re)
if Re<4000
f = 64/Re
else
f = 0.25/(log10(5.74/Re^0.9))^2
end
return f
end
function trapz(t, h)
return 0.5*sum((h[2:end]+h[1:end-1]).*(t[2:end]-t[1:end-1]))
end
function cum_trapz(t, h)
return 0.5*cumsum([0.0; (h[2:end]+h[1:end-1]).*(t[2:end]-t[1:end-1])])
end
function water_exergy(T, p, T0, p0)
if length(T)>1
dh = zeros(length(T))
ds = zeros(length(T))
for i in eachindex(T)
dh[i] = PropsSI("HMASS", "T", T[i], "P", p, "H2O")-
PropsSI("HMASS", "T", T0, "P", p0, "H2O")
ds[i] = PropsSI("SMASS", "T", T[i], "P", p, "H2O")-
PropsSI("SMASS", "T", T0, "P", p0, "H2O")
end
else
dh = PropsSI("HMASS", "T", T, "P", p, "H2O")-
PropsSI("HMASS", "T", T0, "P", p0, "H2O")
ds = PropsSI("SMASS", "T", T, "P", p, "H2O")-
PropsSI("SMASS", "T", T0, "P", p0, "H2O")
end
return dh-T0*ds
end
function water_heat(T, p, T0, p0)
if length(T)>1
dh = zeros(length(T))
ds = zeros(length(T))
for i in eachindex(T)
dh[i] = PropsSI("HMASS", "T", T[i], "P", p, "H2O")-
PropsSI("HMASS", "T", T0, "P", p0, "H2O")
end
else
dh = PropsSI("HMASS", "T", T, "P", p, "H2O")-
PropsSI("HMASS", "T", T0, "P", p0, "H2O")
end
return dh
end
function water_density(T, p)
if length(T)>1
rho = zeros(length(T))
for i in eachindex(T)
rho[i] = PropsSI("D", "T", T[i], "P", p, "H2O")
end
else
rho = PropsSI("D", "T", T, "P", p, "H2O")
end
return rho
end
Out[19]:
In [20]:
PropsSI("D", "T", 300, "P", 1e5, "H2O")
Out[20]:
In [21]:
df_res = CSV.read("results/res_guide.csv")
Out[21]:
In [22]:
CSV.read(df_res[:csv_name][26])
Out[22]:
In [23]:
df_res[:csv_name][17]
Out[23]:
In [24]:
heat_loss = 0.1 # 10% exergy loss
poros = 0.2
# T0 = 20.0+273.15 # [K] reference temperature; defined above
T_in=40.0+273.15 # [K] injection temperature
# p0 = 1.014e5 # [Pa] reference pressure
eta_pump = 0.8; # pump mechanical efficiency
eta_driver = 0.8 # electrical driver efficiency
eta_pp = η_igcc # power plant efficiency
W = 500.0 # [m] width of the reservoir
H = 100.0 # [m] thickness of the reservoir
t_life = 30*365*24*3600.0 # [s] project life time (must be dynamic)
D_well = 0.17 # [m] well diameter
nwells = 2
t_tube = 0.015 # [m] tube thickness
t_cement = 0.050 # [m] cement thickness
ex_dril = 70000.0 # [kJ/m] drilling exergy; practical
depth = 2500.0 # [m]
rho_steel = 7850 # [kg/m^3]
rho_cement = 3150 # [kg/m^3]
ex_steel_th = 6738.0 # 376/0.0558 # [kJ/kg] exergy of Fe
ex_steel = 60000 # [kJ/kg]
ex_cement = 12000 # 3000 # [kJ/kg]
ex_cement_th = 2400 #135/0.056 # [kJ/kg]
ex_methane = ex_gas # 800.0 # [kJ/mol]
cp_water = 4.1855 # [kJ/kg/K]
rho_water = 1000 # [kg/m^3]
mu_water = 1e-3 # [Pa.s]
ex_ccs = 5000 # [kJ/kg CO2]
A_pipe = pi/4*D_well^2 # m^2
c_methane = em_gas # 0.044/ex_methane # [kg CO2/kJ]
Out[24]:
In [25]:
m_steel = nwells*pi/4*((D_well+2*t_tube)^2-D_well^2)*depth*rho_steel # [kg]
m_cement = nwells*pi/4*((D_well+2*t_tube+2*t_cement)^2-(D_well+2*t_tube)^2)*depth*rho_cement # [kg]
Ex_pipe_th=(m_steel*ex_steel_th+m_cement*ex_cement+nwells*depth*ex_dril) # [kJ] theoretical pipe exergy
Ex_pipe_pr=(m_steel*ex_steel+m_cement*ex_cement+nwells*depth*ex_dril) # [kJ] practical pipe exergy
Out[25]:
At which point in time the project starts giving back? Plot the net exergy return vs time, it also shows at which point the project become nonviable. Note: the project life time is the time at which the temperature of the produced water goes below 70 degC. Other criteria can be used for the project life time.
This approach gives a single value of recovery factor for the project life time:
In [26]:
# Just look into one of the simulations
# T_s dp_Pa T_K
fn = "results/Q-50.0-L-650.0-k0.2.csv" # a good one, high perm, low flow
Q = 100/3600 # m^3/s
t1 = CSV.read(fn)
t = float(t1[:T_s]) # [s] simulation time
Th = float(t1[:T_K]) # [K] hot water temperature
# find the project life time
ind_life = findfirst(t1[:T_K].<=(70+273.15))[1]
t_life = t[ind_life]
u_pipe = Q/A_pipe; # [m/s] velocity of water in wells
Re = rho_water*u_pipe*D_well/mu_water # Reynolds number
f = fric_factor(Re)
dp_pipe = nwells*rho_water*f*depth/D_well*u_pipe^2/2 # [Pa] pressure drop in pipes
dp_res = float(t1[:dp_Pa][1:ind_life]) # [Pa] pressure drop in the reservoir
Ex_pump = trapz(t[1:ind_life],Q*(dp_pipe.+dp_res)/1000)
Ex_th = Ex_pump+Ex_pipe_th
Ex_prac = Ex_pump/eta_driver/eta_pp/eta_pump+Ex_pipe_pr # [kJ] practical pumping exergy
heat_prod = Q*rho_water*cp_water*(Th[1:ind_life].-T_in) # [kJ] extracted heat
# Ex_prod = trapz(t,heat_prod.*(1.-T0./Th)) # [kJ] extracted exergy
Ex_prod = Q*trapz(t[1:ind_life], water_density(Th[1:ind_life], p0).*(water_exergy(Th[1:ind_life], p0, T0, p0).-
water_exergy(T_in, p0, T0, p0)))
Ex_prod = Ex_prod/1000 # J to kJ
R_th = (Ex_prod-Ex_th)/Ex_prod # theoretical recovery factor
R_pr = (Ex_prod-Ex_prac)/Ex_prod # theoretical recovery factor
carbon_em = Ex_prac*c_methane # [kg CO2]
Ex_ccs = carbon_em*ex_ccs # [kJ] CCS exergy
R_ze = (Ex_prod-Ex_prac-Ex_ccs)/Ex_prod # zero-emission recovery factor
c = carbon_em/Ex_prod*1000 # [kg CO2/MJ]
Out[26]:
In [27]:
R_th, R_pr, R_ze
Out[27]:
In [28]:
# Just look into one of the simulations
fn = "results/Q-50.0-L-650.0-k0.2.csv" # a good one, high perm, low flow
Q = 100/3600 # m^3/s
t1 = CSV.read(fn)
t = float(t1[:T_s]) # [s] simulation time
Th = float(t1[:T_K]) # [K] hot water temperature
p_out = float(t1[:p_out_Pa])
capital_exergy_th = Ex_pipe_th
capital_exergy_pr = Ex_pipe_pr
# find the project life time
ind_life = length(t) # find(t1[:T_out].<70)[1]
t_life = t[ind_life]
u_pipe = Q/A_pipe; # [m/s] velocity of water in wells
Re = rho_water*u_pipe*D_well/mu_water # Reynolds number
f = fric_factor(Re)
dp_pipe = nwells*rho_water*f*depth/D_well*u_pipe^2/2 # [Pa] pressure drop in pipes
dp_res = float(t1[:dp_Pa][1:ind_life]) # [Pa] pressure drop in the reservoir
dp_dh = water_density(Th[1:ind_life], p0)*9.81*depth-p_out[1:ind_life]
Ex_pump = cum_trapz(t[1:ind_life],Q*((dp_pipe.+dp_res)+dp_dh)/1000)
Ex_th = Ex_pump.+Ex_pipe_th
Ex_prac = Ex_pump/eta_driver/eta_pp/eta_pump.+Ex_pipe_pr # [kJ] practical pumping exergy
heat_prod = Q*rho_water*cp_water*(Th[1:ind_life].-T_in) # [kJ] extracted heat
# Ex_prod = trapz(t,heat_prod.*(1-T0./Th)) # [kJ] extracted exergy
Ex_prod = Q*cum_trapz(t[1:ind_life], water_density(Th[1:ind_life], p0).*(water_exergy(Th[1:ind_life], p0, T0, p0).-
water_exergy(T_in, p0, T0, p0)))
Ex_prod = Ex_prod/1000 # J to kJ
R_th = (Ex_prod-Ex_th)./Ex_prod # theoretical recovery factor
R_pr = (Ex_prod-Ex_prac)./Ex_prod # theoretical recovery factor
carbon_em = Ex_prac*c_methane # [kg CO2]
Ex_ccs = carbon_em*ex_ccs # [kJ] CCS exergy
R_ze = (Ex_prod-Ex_prac-Ex_ccs)./Ex_prod # zero-emission recovery factor
c = carbon_em/Ex_prod*1000 # [kg CO2/MJ]
Out[28]:
In [30]:
figure(figsize=(4,3))
plot(float(t[1:ind_life])/(365*24*3600),R_th, "--", label="Theoretical")
plot(float(t[1:ind_life])/(365*24*3600),R_pr, label="Practical")
plot(float(t[1:ind_life])/(365*24*3600),R_ze, "-.", label="Zero-emission")
plot([0, 30], [0, 0], "k--")
legend()
# grid()
xlabel("Production time [year]")
ylabel("Recovery factor")
#yscale("log")
axis([0,30, -1,1])
tight_layout()
savefig("figs/rec_fac_time.png")
savefig("figs/rec_fac_time.svg")
# cumsum(rand(10))
The point at which the recovery factor is zero shows the time at which all the initial exergy consumption is returned (or paid back) by the geothermal reservoir.
What was missing from my previous calculation is the down hole pump for the extraction of water (note that the bottomhole pressure is less than the hydrostatic pressure). The pump must overcome the difference between the hydrostatic pressure and the bottomhole pressure, i.e., $$ W_{pump,dh} = Q_{water}(\rho_{water}gh_{well}-p_{out})$$ where $p_{out}$ is the production well downhole pressure. This needs to be included in the analysis.
In [36]:
df_res[:csv_name]
Out[36]:
In [41]:
# create empty data frames to save the analyzed files
df_geothermal = DataFrame(q_m3_s=[], k_m2=[], L_m=[], t_life_s=[], R_th=[], R_pr=[],
R_ze=[], c_CO2=[], COP=[], c_CO2_heat=[], heat_per_m3=[])
i = 0
for fn in df_res[:csv_name]
i=i+1
Q = df_res[:q_m3_h][i]/3600 # m3/s
perm = df_res[:k_D][i]*1e-12 # m2
L = df_res[:L_m][i] # m
println(fn)
t1 = CSV.read(fn)
# deleterows!(t1,find(isna(t1[:,:Time_s_])))
ind_life = findfirst(t1[:T_K].<=(70+273.15))[1]
t=float(t1[:T_s][1:ind_life]) # [s] simulation time
t_life = t[ind_life]
Th = float(t1[:T_K][1:ind_life]) # [K] hot water temperature
p_out = float(t1[:p_out_Pa][1:ind_life])
dp_res = float(t1[:dp_Pa][1:ind_life]) # [Pa] pressure drop in the reservoir
u_pipe = Q/A_pipe # [m/s] velocity of water in wells
Re = rho_water*u_pipe*D_well/mu_water # Reynolds number
f=fric_factor(Re)
dp_pipe = nwells*rho_water*f*depth/D_well*u_pipe^2/2 # [Pa] pressure drop in pipes
dp_dh = water_density(Th[1:ind_life], p0)*9.81*depth-p_out # [Pa]
Ex_pump = trapz(t,Q*(dp_pipe.+dp_res.+dp_dh)/1000)
Ex_th = Ex_pump+Ex_pipe_th
Ex_prac = Ex_pump/eta_driver/eta_pp/eta_pump+Ex_pipe_pr # [kJ] practical pumping exergy
# heat_prod = Q*rho_water*cp_water*(Th-T_in) # [kJ] extracted heat
# Ex_prod = trapz(t,heat_prod.*(1-T0./Th)) # [kJ] extracted exergy
Ex_prod = Q*trapz(t, water_density(Th, p0).*(water_exergy(Th, p0, T0, p0).-
water_exergy(T_in, p0, T0, p0)))
heat_prod = Q*trapz(t, water_density(Th, p0).*(water_heat(Th, p0, T0, p0).-
water_heat(T_in, p0, T0, p0)))
heat_prod = heat_prod/1000 # J to kJ
heat_waste = heat_prod*heat_loss
Ex_prod = Ex_prod/1000 # J to kJ
Ex_waste = Ex_prod*heat_loss
R_th = (Ex_prod-Ex_waste-Ex_th)/Ex_prod # theoretical recovery factor
R_pr = (Ex_prod-Ex_waste-Ex_prac)/Ex_prod # theoretical recovery factor
carbon_em = Ex_prac*c_methane # [kg CO2]
Ex_ccs = carbon_em*ex_ccs # [kJ] CCS exergy
R_ze = (Ex_prod-Ex_waste-Ex_prac-Ex_ccs)/Ex_prod # zero-emission recovery factor
c = carbon_em/(Ex_prod-Ex_waste)*1000 # [kg CO2/MJ]
c_heat = carbon_em/(heat_prod-heat_waste)*1000 # [kg CO2/MJ]
cop = (heat_prod-heat_waste)/Ex_prac # kJ/kJ
push!(df_geothermal, [Q, perm, L, t_life, R_th, R_pr, R_ze, c, cop, c_heat, (heat_prod-heat_loss)/(Q*t_life)])
end
In [46]:
df_geothermal
Out[46]:
In [47]:
CSV.write("geothermal_final_results.csv", df_geothermal)
Out[47]:
In [32]:
(df_geothermal)
Out[32]:
In [44]:
plot(df_geothermal[:L_m], df_geothermal[:t_life_s]/365/24/3600, "o")
Out[44]:
In [42]:
f_name = ["Geothermal_Q150_2Q150-20mD-L800.csv"
"Geothermal_Q150_2Q150-50mD-L800.csv"
"Geothermal_Q150_2Q150-100mD-L800.csv"
"Geothermal_Q150_2Q150-500mD-L800.csv"
"Geothermal_Q150_2Q150-1D-L800.csv"]
i=0
R_th=zeros(5)
R_pr=zeros(5)
R_ze=zeros(5)
c=zeros(5)
for fn in f_name
i=i+1
Q=100/3600 # m^3/s
t1 = CSV.read(fn)
t=float(t1[:Time_s_]) # [s] simulation time
Th = float(t1[:T_out])+273.15 # [K] hot water temperature
u_pipe = Q/A_pipe; # [m/s] velocity of water in wells
Re = rho_water*u_pipe*D_well/mu_water # Reynolds number
f=fric_factor(Re)
dp_pipe = nwells*rho_water*f*depth/D_well*u_pipe^2/2 # [Pa] pressure drop in pipes
dp_res = float(t1[:Dp_Mpa_])*1e6 # [Pa] pressure drop in the reservoir
Ex_pump = trapz(t,Q*(dp_pipe+dp_res)/1000)
Ex_th = Ex_pump+Ex_pipe_th
Ex_prac = Ex_pump/eta_driver/eta_pp/eta_pump+Ex_pipe_pr # [kJ] practical pumping exergy
heat_prod = Q*rho_water*cp_water*(Th-T_in) # [kJ] extracted heat
# Ex_prod = trapz(t,heat_prod.*(1-T0./Th)) # [kJ] extracted exergy
Ex_prod = Q*rho_water*trapz(t, water_exergy(Th, p0, T0, p0)-
water_exergy(T_in, p0, T0, p0))
Ex_prod = Ex_prod/1000 # J to kJ
R_th[i] = (Ex_prod-Ex_th)/Ex_prod # theoretical recovery factor
R_pr[i] = (Ex_prod-Ex_prac)/Ex_prod # theoretical recovery factor
carbon_em = Ex_prac*c_methane # [kg CO2]
Ex_ccs = carbon_em*ex_ccs # [kJ] CCS exergy
R_ze[i] = (Ex_prod-Ex_prac-Ex_ccs)/Ex_prod # zero-emission recovery factor
c[i] = carbon_em/Ex_prod*1000 # [kg CO2/MJ]
end
In [43]:
figure(figsize=(5,3))
perm
plot(perm_range, R_th, "-o", perm_range, R_pr, "--^k", perm_range, R_ze, "-vr")
xlabel("Permeability [mD]", fontsize=12)
ylabel("Recovery factor [-]", fontsize=12)
legend(["Theoretical", "Practical", "Zero-emission"], loc=4, fontsize=12)
plot([0, 1000], [0, 0], ":k")
# text(800, 0.25, "(a)", fontsize=14)
axis([0, 1000, -1.0, 1.0])
savefig("perm_effect.png")
In [49]:
figure(figsize=(5,3))
plot(perm_range, c, "-o")
xlabel("Permeability [mD]", fontsize=12)
ylabel("CO2 emission [kg CO2/MJ]", fontsize=12)
#legend(["Theoretical", "Practical", "Zero-emission"], loc=4, fontsize=12)
axis([0, 1000, 0.0, 0.3])
savefig("perm_effect_CO2_emission.svg")
In [31]:
f_name = ["Geothermal_Q250_Q250-20mD-L800.csv"
"Geothermal_Q250_Q250-50mD-L800.csv"
"Geothermal_Q250_Q250-100mD-L800.csv"
"Geothermal_Q250_Q250-500mD-L800.csv"
"Geothermal_Q250_Q250-1D-L800.csv"]
i=0
R_th=zeros(5)
R_pr=zeros(5)
R_ze=zeros(5)
c=zeros(5)
for fn in f_name
i=i+1
Q=250/3600 # m^3/s
t1 = CSV.read(fn)
t=float(t1[:Time_s_]) # [s] simulation time
Th = float(t1[:T_out])+273.15 # [K] hot water temperature
u_pipe = Q/A_pipe; # [m/s] velocity of water in wells
Re = rho_water*u_pipe*D_well/mu_water # Reynolds number
f=fric_factor(Re)
dp_pipe = nwells*rho_water*f*depth/D_well*u_pipe^2/2 # [Pa] pressure drop in pipes
dp_res = float(t1[:Dp_Mpa_])*1e6 # [Pa] pressure drop in the reservoir
Ex_pump = trapz(t,Q*(dp_pipe+dp_res)/1000)
Ex_th = Ex_pump+Ex_pipe_th
Ex_prac = Ex_pump/eta_driver/eta_pp/eta_pump+Ex_pipe_pr # [kJ] practical pumping exergy
heat_prod = Q*rho_water*cp_water*(Th-T_in) # [kJ] extracted heat
Ex_prod = Q*rho_water*trapz(t, water_exergy(Th, p0, T0, p0)-
water_exergy(T_in, p0, T0, p0))
Ex_prod = Ex_prod/1000 # J to kJ
# Ex_prod = trapz(t,heat_prod.*(1-T0./Th)) # [kJ] extracted exergy
R_th[i] = (Ex_prod-Ex_th)/Ex_prod # theoretical recovery factor
R_pr[i] = (Ex_prod-Ex_prac)/Ex_prod # theoretical recovery factor
carbon_em = Ex_prac*c_methane # [kg CO2]
Ex_ccs = carbon_em*ex_ccs # [kJ] CCS exergy
R_ze[i] = (Ex_prod-Ex_prac-Ex_ccs)/Ex_prod # zero-emission recovery factor
c[i] = carbon_em/Ex_prod*1000 # [kg CO2/MJ]
end
In [32]:
figure(figsize=(5,3))
plot(perm_range, R_th, "-o", perm_range, R_pr, "--^k", perm_range, R_ze, "-vr")
xlabel("Permeability [mD]", fontsize=12)
ylabel("Recovery factor [-]", fontsize=12)
legend(["Theoretical", "Practical", "Zero-emission"], loc=4, fontsize=12)
plot([0, 1000], [0, 0], ":k")
text(800, 0.25, "(b)", fontsize=14)
axis([0, 1000, -1.0, 1.0])
savefig("perm_effect_high_flow.svg")
In [33]:
figure(figsize=(5,3))
plot(perm_range, c, "-o")
xlabel("Permeability [mD]", fontsize=14)
ylabel("CO2 emission [kg CO2/MJ]", fontsize=14)
#legend(["Theoretical", "Practical", "Zero-emission"], loc=4, fontsize=12)
savefig("perm_effect_CO2_emission_high_flow.svg")
In [136]:
f_name = ["Geothermal_Q100_Q100-100mD-L470.csv", "Geothermal_Q150_Q150-100mD-L600.csv",
"Geothermal_Q200_Q200-100mD-L700.csv", "Geothermal_Q250_Q250-100mD-L800.csv"]
i=0
R_th=zeros(4)
R_pr=zeros(4)
R_ze=zeros(4)
c=zeros(4)
for fn in f_name
i=i+1
Q=Q_range[i] # m^3/s
t1 = CSV.read(fn)
t=float(t1[:Time_s_]) # [s] simulation time
Th = float(t1[:T_out])+273.15 # [K] hot water temperature
u_pipe = Q/A_pipe; # [m/s] velocity of water in wells
Re = rho_water*u_pipe*D_well/mu_water # Reynolds number
f=fric_factor(Re)
dp_pipe = nwells*rho_water*f*depth/D_well*u_pipe^2/2 # [Pa] pressure drop in pipes
dp_res = float(t1[:Dp_Mpa_])*1e6 # [Pa] pressure drop in the reservoir
Ex_pump = trapz(t,Q*(dp_pipe+dp_res)/1000)
Ex_th = Ex_pump+Ex_pipe_th
Ex_prac = Ex_pump/eta_driver/eta_pp/eta_pump+Ex_pipe_pr # [kJ] practical pumping exergy
heat_prod = Q*rho_water*cp_water*(Th-T_in) # [kJ] extracted heat
Ex_prod = Q*rho_water*trapz(t, water_exergy(Th, p0, T0, p0)-
water_exergy(T_in, p0, T0, p0))
Ex_prod = Ex_prod/1000 # J to kJ
#Ex_prod = trapz(t,heat_prod.*(1-T0./Th)) # [kJ] extracted exergy
R_th[i] = (Ex_prod-Ex_th)/Ex_prod # theoretical recovery factor
R_pr[i] = (Ex_prod-Ex_prac)/Ex_prod # theoretical recovery factor
carbon_em = Ex_prac*c_methane # [kg CO2]
Ex_ccs = carbon_em*ex_ccs # [kJ] CCS exergy
R_ze[i] = (Ex_prod-Ex_prac-Ex_ccs)/Ex_prod # zero-emission recovery factor
c[i] = carbon_em/Ex_prod*1000 # [kg CO2/MJ]
end
In [137]:
figure(figsize=(5,3))
plot(L_range, Q_range*3600, "-o")
xlabel("Well spacing [m]", fontsize=12)
ylabel("Flow rate "*L"[m^3/h]", fontsize=12)
savefig("well_space_flow_rate.svg")
In [138]:
figure(figsize=(5,3))
plot(L_range, R_th, "-o", L_range, R_pr, "--^k", L_range, R_ze, "-vr")
xlabel("Well spacing [m]", fontsize=12)
ylabel("Recovery factor [-]", fontsize=12)
legend(["Theoretical", "Practical", "Zero-emission"], loc=3, fontsize=12)
savefig("well_space.svg")
In [139]:
f_name = ["Geothermal_Q100_Q100-500mD-L470.csv", "Geothermal_Q150_Q150-500mD-L600.csv",
"Geothermal_Q200_Q200-500mD-L700.csv", "Geothermal_Q250_Q250-500mD-L800.csv"]
i=0
R_th=zeros(4)
R_pr=zeros(4)
R_ze=zeros(4)
c=zeros(4)
for fn in f_name
i=i+1
Q=Q_range[i] # m^3/s
t1 = CSV.read(fn)
t=float(t1[:Time_s_]) # [s] simulation time
Th = float(t1[:T_out])+273.15 # [K] hot water temperature
u_pipe = Q/A_pipe; # [m/s] velocity of water in wells
Re = rho_water*u_pipe*D_well/mu_water # Reynolds number
f=fric_factor(Re)
dp_pipe = nwells*rho_water*f*depth/D_well*u_pipe^2/2 # [Pa] pressure drop in pipes
dp_res = float(t1[:Dp_Mpa_])*1e6 # [Pa] pressure drop in the reservoir
Ex_pump = trapz(t,Q*(dp_pipe+dp_res)/1000)
Ex_th = Ex_pump+Ex_pipe_th
Ex_prac = Ex_pump/eta_driver/eta_pp/eta_pump+Ex_pipe_pr # [kJ] practical pumping exergy
heat_prod = Q*rho_water*cp_water*(Th-T_in) # [kJ] extracted heat
# Ex_prod = trapz(t,heat_prod.*(1-T0./Th)) # [kJ] extracted exergy
Ex_prod = Q*rho_water*trapz(t, water_exergy(Th, p0, T0, p0)-
water_exergy(T_in, p0, T0, p0))
Ex_prod = Ex_prod/1000 # J to kJ
R_th[i] = (Ex_prod-Ex_th)/Ex_prod # theoretical recovery factor
R_pr[i] = (Ex_prod-Ex_prac)/Ex_prod # theoretical recovery factor
carbon_em = Ex_prac*c_methane # [kg CO2]
Ex_ccs = carbon_em*ex_ccs # [kJ] CCS exergy
R_ze[i] = (Ex_prod-Ex_prac-Ex_ccs)/Ex_prod # zero-emission recovery factor
c[i] = carbon_em/Ex_prod*1000 # [kg CO2/MJ]
end
In [140]:
figure(figsize=(5,3))
plot(L_range, R_th, "-o", L_range, R_pr, "--^k", L_range, R_ze, "-vr")
xlabel("Well spacing [m]", fontsize=14)
ylabel("Recovery factor [-]", fontsize=14)
legend(["Theoretical", "Practical", "Zero-emission"], loc=3, fontsize=12)
savefig("well_space_500.svg")
In [141]:
f_name = ["Geothermal_Q100_Q100-50mD-L470.csv", "Geothermal_Q150_Q150-50mD-L600.csv",
"Geothermal_Q200_Q200-50mD-L700.csv", "Geothermal_Q250_Q250-50mD-L800.csv"]
i=0
R_th=zeros(4)
R_pr=zeros(4)
R_ze=zeros(4)
c=zeros(4)
for fn in f_name
i=i+1
Q=Q_range[i] # m^3/s
t1 = CSV.read(fn)
t=float(t1[:Time_s_]) # [s] simulation time
Th = float(t1[:T_out])+273.15 # [K] hot water temperature
u_pipe = Q/A_pipe; # [m/s] velocity of water in wells
Re = rho_water*u_pipe*D_well/mu_water # Reynolds number
f=fric_factor(Re)
dp_pipe = nwells*rho_water*f*depth/D_well*u_pipe^2/2 # [Pa] pressure drop in pipes
dp_res = float(t1[:Dp_Mpa_])*1e6 # [Pa] pressure drop in the reservoir
Ex_pump = trapz(t,Q*(dp_pipe+dp_res)/1000)
Ex_th = Ex_pump+Ex_pipe_th
Ex_prac = Ex_pump/eta_driver/eta_pp/eta_pump+Ex_pipe_pr # [kJ] practical pumping exergy
heat_prod = Q*rho_water*cp_water*(Th-T_in) # [kJ] extracted heat
# Ex_prod = trapz(t,heat_prod.*(1-T0./Th)) # [kJ] extracted exergy
Ex_prod = Q*rho_water*trapz(t, water_exergy(Th, p0, T0, p0)-
water_exergy(T_in, p0, T0, p0))
Ex_prod = Ex_prod/1000 # J to kJ
R_th[i] = (Ex_prod-Ex_th)/Ex_prod # theoretical recovery factor
R_pr[i] = (Ex_prod-Ex_prac)/Ex_prod # theoretical recovery factor
carbon_em = Ex_prac*c_methane # [kg CO2]
Ex_ccs = carbon_em*ex_ccs # [kJ] CCS exergy
R_ze[i] = (Ex_prod-Ex_prac-Ex_ccs)/Ex_prod # zero-emission recovery factor
c[i] = carbon_em/Ex_prod*1000 # [kg CO2/MJ]
end
In [142]:
figure(figsize=(5,3))
plot(L_range, R_th, "-o", L_range, R_pr, "--^k", L_range, R_ze, "-vr")
xlabel("Well spacing [m]", fontsize=14)
ylabel("Recovery factor [-]", fontsize=14)
legend(["Theoretical", "Practical", "Zero-emission"], loc=3, fontsize=12)
savefig("well_space_50.svg")
m
i=0;
%for perm = 3e-14:0.1e-14:10e-14
for Q = 100/3600:100/3600:1000/3600
i=i+1;
% calculation starts here ----------------------------------------------->
m_steel = nwells*pi()/4*((D_well+2*t_tube)^2-D_well^2)*depth*rho_steel; % [kg]
m_cement = nwells*pi()/4* ...
((D_well+2*t_tube+2*t_cement)^2-(D_well+2*t_tube)^2)*depth*rho_cement; % [kg]
A_res = W*H; % [m^2] reservoir cross section
A_pipe = pi()/4*D_well^2;
u_res = Q/A_res; % [m/s] Darcy velocity
u_pipe = Q/A_pipe; % [m/s] velocity of water in wells
Re = rho_water*u_pipe*D_well/mu_water; % Reynolds number
dp_pipe = nwells*rho_water*f*depth/D_well*u_pipe^2/2; % [Pa] pressure drop in pipes
dp_res = mu_water/perm*u_res*L; % [Pa] pressure drop in the reservoir
Ex_pump = Q*(dp_pipe+dp_res)/1000;
Ex_pipe_th=(m_steel*ex_steel_th+m_cement*ex_cement+ ...
nwells*depth*ex_dril)/t_life; % [kJ] theoretical pipe exergy
Ex_th = Ex_pump+Ex_pipe_th;
Ex_pipe_pr=(m_steel*ex_steel+m_cement*ex_cement+ ...
nwells*depth*ex_dril)/t_life; % [kJ] practical pipe exergy
Ex_prac = Ex_pump/eta_driver/eta_pp/eta_pump+Ex_pipe_pr; % [kJ] practical pumping exergy
heat_prod = Q*rho_water*cp_water*(Th-T0) # [kJ] extracted heat
Ex_prod = trapz(t,heat_prod.*(1-T0./Th)) # [kJ] extracted exergy
R_th(i) = (Ex_prod-Ex_th)/Ex_prod; % theoretical recovery factor
R_pr(i) = (Ex_prod-Ex_prac)/Ex_prod; % theoretical recovery factor
carbon_em = Ex_prac*c_methane; % [kg CO2]
Ex_ccs = carbon_em*ex_ccs; % [kJ] CCS exergy
R_ze(i) = (Ex_prod-Ex_prac-Ex_ccs)/Ex_prod; % zero-emission recovery factor
c(i) = carbon_em/Ex_prod*1000; % [kg CO2/MJ]
% end of calculation ----------------------------------------------->
end
%plot(perm, R_th, 'o', perm, R_pr, '^');
Q = 100/3600:100/3600:1000/3600;
plot(Q, R_th, Q, R_th, 'o', ...
Q, R_pr, Q, R_pr, '^', ...
Q, R_ze, Q, R_ze, 's');
xlabel('flow rate (m^3/s)');
perm = 3e-14:0.1e-14:10e-14;
%plot(perm, R_th, perm, R_th, 'o', ...
% perm, R_pr, perm, R_pr, '^', ...
% perm, R_ze, perm, R_ze, 's');
%xlabel('Permeability [m^2]');
ylabel('Recovery factors');
legend('Theoretical', '','Practical', '', 'Zero-emission','');
In [ ]:
In [87]:
# file_list = readdir()
# f_list = Array{Any}(20)
# j=0
# for i in eachindex(file_list)
# if file_list[i][end-2:end]=="csv"
# j+=1
# f_list[j] = file_list[i]
# end
# end
f_list =
["Geothermal_Q100_Q100-100mD-L470.csv"
"Geothermal_Q100_Q100-1D-L470.csv"
"Geothermal_Q100_Q100-20mD-L470.csv"
"Geothermal_Q100_Q100-500mD-L470.csv"
"Geothermal_Q100_Q100-50mD-L470.csv"
"Geothermal_Q150_2Q150-100mD-L800.csv"
"Geothermal_Q150_2Q150-1D-L800.csv"
"Geothermal_Q150_2Q150-20mD-L800.csv"
"Geothermal_Q150_2Q150-500mD-L800.csv"
"Geothermal_Q150_2Q150-50mD-L800.csv"
"Geothermal_Q150_Q150-100mD-L600.csv"
"Geothermal_Q150_Q150-1D-L600.csv"
"Geothermal_Q150_Q150-20mD-L600.csv"
"Geothermal_Q150_Q150-500mD-L600.csv"
"Geothermal_Q150_Q150-50mD-L600.csv"
"Geothermal_Q200_Q200-100mD-L700.csv"
"Geothermal_Q200_Q200-1D-L700.csv"
"Geothermal_Q200_Q200-20mD-L700.csv"
"Geothermal_Q200_Q200-500mD-L700.csv"
"Geothermal_Q200_Q200-50mD-L700.csv"
"Geothermal_Q250_Q250-100mD-L800.csv"
"Geothermal_Q250_Q250-1D-L800.csv"
"Geothermal_Q250_Q250-20mD-L800.csv"
"Geothermal_Q250_Q250-500mD-L800.csv"
"Geothermal_Q250_Q250-50mD-L800.csv"]
# flow[m3/h] perm[mD] L[m]
sim_params =
[100 100 470
100 1000 470
100 20 470
100 500 470
100 50 470
150 100 600
150 1000 600
150 20 600
150 500 600
150 50 600
150 100 800
150 1000 800
150 20 800
150 500 800
150 50 800
200 100 700
200 1000 700
200 20 700
200 500 700
200 50 700
250 100 800
250 1000 800
250 20 800
250 500 800
250 50 800]
Out[87]:
In [26]:
find(sim_params[:,2].==100)
length(f_list)
Out[26]:
In [92]:
i = 0
n_data = length(f_list)
R_th = zeros(n_data)
R_pr = zeros(n_data)
R_ze = zeros(n_data)
Q = sim_params[:, 1]/3600 # m3/s
perm = sim_params[:, 2]/1000*1e-12 # m2
L = sim_params[:, 3] # m
c = zeros(n_data)
for fn in f_list
i=i+1
println(fn)
t1 = CSV.read(fn)
deleterows!(t1,find(isna(t1[:,:Time_s_])))
t=float(t1[:Time_s_]) # [s] simulation time
Th = float(t1[:T_out])+273.15 # [K] hot water temperature
u_pipe = Q[i]/A_pipe; # [m/s] velocity of water in wells
Re = rho_water*u_pipe*D_well/mu_water # Reynolds number
f=fric_factor(Re)
dp_pipe = nwells*rho_water*f*depth/D_well*u_pipe^2/2 # [Pa] pressure drop in pipes
dp_res = float(t1[:Dp_Mpa_])*1e6 # [Pa] pressure drop in the reservoir
Ex_pump = trapz(t,Q[i]*(dp_pipe+dp_res)/1000)
Ex_th = Ex_pump+Ex_pipe_th
Ex_prac = Ex_pump/eta_driver/eta_pp/eta_pump+Ex_pipe_pr # [kJ] practical pumping exergy
heat_prod = Q[i]*rho_water*cp_water*(Th-T_in) # [kJ] extracted heat
# Ex_prod = trapz(t,heat_prod.*(1-T0./Th)) # [kJ] extracted exergy
Ex_prod = Q[i]*rho_water*trapz(t, water_exergy(Th, p0, T0, p0)-
water_exergy(T_in, p0, T0, p0))
Ex_prod = Ex_prod/1000 # J to kJ
R_th[i] = (Ex_prod-Ex_th)/Ex_prod # theoretical recovery factor
R_pr[i] = (Ex_prod-Ex_prac)/Ex_prod # theoretical recovery factor
carbon_em = Ex_prac*c_methane # [kg CO2]
Ex_ccs = carbon_em*ex_ccs # [kJ] CCS exergy
R_ze[i] = (Ex_prod-Ex_prac-Ex_ccs)/Ex_prod # zero-emission recovery factor
c[i] = carbon_em/Ex_prod*1000 # [kg CO2/MJ]
end
In [93]:
# Q=250, L = 800
ind = [11:15; 21:25]
spl_R_pr = Spline2D(perm[ind], Q[ind], R_pr[ind], kx=1, ky=1, s=16)
Qi = collect(linspace(minimum(Q[ind]), maximum(Q[ind]), 20))
permi = collect(linspace(minimum(perm[ind]), maximum(perm[ind]), 20))
R_pri = evalgrid(spl_R_pr, permi, Qi)
Out[93]:
In [94]:
plot3D(perm[ind], Q[ind], R_pr[ind], "o")
surf(permi, Qi, R_pri)
xlabel("Permeability [m^2]")
ylabel("flow rate [m^3/s]")
zlabel("Recovery factor (pr)")
Out[94]: