In [4]:
using Dierckx, DataFrames, PyPlot, JLD, CoolProp, CSV
include("SetPyPlot.jl")
SetPyPlot.setrcparam()
@load "input_param.jld"
Out[4]:
This data is obtained from the following webpage: Source: http://nlog.nl/reservoireigenschappen
In [29]:
phi_data = CSV.read("TNO_data.csv")
perm_data = dropmissing(phi_data, :Perm_arith)
Out[29]:
In [35]:
figure(figsize=(4,3))
plot(phi_data.Basis_tvd_m, phi_data.Porositeit, "o")
xlabel("Depth [m]")
ylabel("Porosity [vol%]")
Out[35]:
In [34]:
# plotyy(x1, y1, x2, y2; fig_size, x_label, y1_label, y2_label, y1_style, y2_style)
SetPyPlot.plotyy(phi_data.Basis_tvd_m, phi_data.Porositeit,
perm_data.Basis_tvd_m, perm_data.Perm_arith, y1_style = "ob", y2_style="vr",
x_label="Depth [m]", y1_label="Porosity [vol%]", y2_label="mili Darcy")
# plot(phi_data.Basis_tvd_m, phi_data.Perm_arith)
Out[34]:
In [13]:
em_gas*1e3
Out[13]:
In [3]:
df_geo = CSV.read("geothermal_final_results.csv")
Out[3]:
In [4]:
flow_rates = df_geo[:q_m3_s] # [m^3/s]
k = df_geo[:k_m2] # [m^2]
L = df_geo[:L_m] # [m]
flow_range = union(flow_rates)
k_range = union(k)
L_range = union(L)
Out[4]:
In [5]:
k_range
Out[5]:
In [6]:
line_style = ["-", "--", "-.", "-", ":"]
for L_well in L_range
figure(figsize=(5,3))
xlabel("Permeability [m\$^2\$]")
ylabel("Practical recovery factor [-]")
i=0
for q in flow_range
i+=1
df_temp = df_geo[(df_geo[:q_m3_s].==q) .& (df_geo[:L_m].==L_well), :]
plot(df_temp[:k_m2], df_temp[:R_pr], line_style[i], label = "q = $(q*3600) m\$^3\$/h")
end
legend()
plot([0, 1e-12], [0,0], "--k")
axis([0, 1e-12, -1,1])
# title("L = $L_well [m]")
savefig("figs/practical_recovery_perm.png")
end
In [7]:
for L_well in L_range
figure(figsize=(4,3))
xlabel("Permeability [m\$^2\$]")
ylabel("Zero-emission recovery factor [-]")
for q in flow_range
df_temp = df_geo[(df_geo[:q_m3_s].==q) .& (df_geo[:L_m].==L_well), :]
plot(df_temp[:k_m2], df_temp[:R_ze], line_style[i], label = "q = $(q*3600) m\$^3\$/h")
# plot(df_temp[:k_m2], df_temp[:R_pr], label = "q = $(q*3600) m\$^3\$/h")
end
legend()
plot([0, 1e-12], [0,0], "--k")
axis([0, 1e-12, -1,1])
title("L = $L_well [m]")
end
In [8]:
for k_res in k_range
figure(figsize=(5,3))
xlabel("Flow rate \L[m^3/h]")
ylabel("Practical recovery factor [-]")
for L_well in L_range
df_temp = df_geo[(df_geo[:k_m2].==k_res) & (df_geo[:L_m].==L_well), :]
plot(df_temp[:q_m3_s]*3600, df_temp[:R_pr], line_style[i], alpha = 0.5, label = "L = $L_well [m]")
end
legend()
#plot([0, 1e-12], [0,0], "--k")
#axis([0, 1e-12, -1,1])
title("k = $(k_res*1e12) [D]")
end
In [9]:
for k_res in k_range
figure(figsize=(4,3))
xlabel("Flow rate [m\$^3\$/h]")
ylabel("Project life time [year]")
for L_well in L_range
df_temp = df_geo[(df_geo[:k_m2].==k_res) .& (df_geo[:L_m].==L_well), :]
plot(df_temp[:q_m3_s]*3600, df_temp[:t_life_s]/(3600*24*365), line_style[i], alpha = 0.8, label = "L = $L_well [m]")
end
legend()
#plot([0, 1e-12], [0,0], "--k")
#axis([0, 1e-12, -1,1])
# plot(df_geo[:q_m3_s]*3600, df_geo[:t_life_s]/(3600*24*365), "o")
savefig("figs/project_life_time.png")
# title("k = $(k_res*1e12) [D]")
end
In [10]:
figure(figsize=(4,3))
semilogy(df_geo[:q_m3_s]*3600, df_geo[:c_CO2], alpha = 0.5, "o")
semilogy([0, 260], [em_gas, em_gas]*1e3, "--k")
xlabel("Flow rate [m\$^3\$/h]")
ylabel("CO\$_2\$ emission [kg CO\$_2\$/MJ]")
axis([0, 260, 1e-2, 30])
savefig("figs/CO2_emission_flow.png")
In [11]:
figure(figsize=(4,3))
semilogy(df_geo[:k_m2], df_geo[:c_CO2], alpha = 0.5, "o")
semilogy([0, 1.1e-12], [em_gas, em_gas]*1e3, "--k")
xlabel("Permeability [m\$^2\$]")
ylabel("CO\$_2\$ emission [kg CO\$_2\$/MJ]")
axis([0, 1.1e-12, 1e-2, 30])
savefig("figs/CO2_emission_perm.png")
In [26]:
figure(figsize=(4,3))
semilogy(df_geo[:k_m2], df_geo[:c_CO2]*0.15, alpha = 0.5, "o")
semilogy([0, 1.1e-12], [em_gas, em_gas]*1e3, "--k")
xlabel("Permeability [m\$^2\$]")
ylabel("CO\$_2\$ emission [kg CO\$_2\$/MJ]")
savefig("figs/CO2_emission_perm.png")
In [12]:
df_geo.colindex
In [13]:
T_c = 2.0+273.15 # [K]
T_h = 22+273.15 # [K]
COP_th = T_h/(T_h-T_c)
COP_air = 2.0 # average value (per unit electricity)
COP_geo = 5.0 # shallow heat (comes from the sun)
COP_methane = 1.0 # ignoring heat loss
t_life_crit = 30 # years
c_CO2_crit = em_gas*1e3*2.0 # kg/MJ
u_ht = 1.0 # W/m2/K
A_ht = 4*10.0*3.0 # only surrounding walls
Q_heat = u_ht*A_ht*(T_h-T_c)/1000 # kJ/s
n_houses = 400 # number of houses to be covered by the geothermal project
m_water = Q_heat*1000*n_houses/(PropsSI("HMASS", "T", 80+273.15, "P", 1e5, "H2O")-PropsSI("HMASS", "T", 40+273.15, "P", 1e5, "H2O"))
q_water = m_water/PropsSI("D", "T", 40+273.15, "P", 1e5, "H2O") # [m3/s]
df_practical = df_geo[(df_geo[:t_life_s].>(30.0*24*3600)) & (df_geo[:R_ze].>0.0) & (df_geo[:c_CO2].<c_CO2_crit), :]
In [14]:
flow_rates_pr = df_practical[:q_m3_s] # [m^3/s]
k_pr = df_practical[:k_m2] # [m^2]
L_pr = df_practical[:L_m] # [m]
flow_range_pr = union(flow_rates_pr)
k_range_pr = union(k_pr)
L_range_pr = union(L_pr)
In [15]:
k_range_pr
In [16]:
figure(figsize=(5,3))
plot(df_practical[:q_m3_s]*3600, df_practical[:c_CO2], alpha = 0.5, "o")
plot([0, 260], [em_gas, em_gas]*1e3, "--k")
xlabel("Flow rate [m^3/h]")
ylabel("CO2 emission [kg CO2/MJ]")
tight_layout()
savefig("carbon_emission_geothermal.png")
The theoretical COP can be calculated based on the temperature of the outside (temperature outside of the house, where heat is extracted) and the desirable temperature inside (where heat is delivered). It is calculated by $$COP_{th} = \frac{T_h}{T_h-T_c}$$
In [17]:
df_current = df_geo[df_geo[:COP].<1.0, :]
Out[17]:
:MW_CO2
:MW_coal
:MW_gas
:e_igcc
:e_ngcc
:e_pc
:em_coal
:em_gas
:eta_driver
:eta_pump
:ex_ccs
:ex_coal
:ex_gas
:η_igcc
:η_ngcc
:η_pc
In [18]:
figure(figsize=(5,3))
plot(df_geo[:k_m2]*1e12, df_geo[:COP], "o", alpha = 0.5)
plot([0,1.1], [1.0, 1.0], "--k", label = "methane")
plot([0,1.1], [2.0, 2.0], "-r", label = "heat pump")
plot([0,1.1], [5.0, 5.0], ":g", label = "geo heat pump")
axis([0,1.05,0,40])
legend()
ylabel("COP")
xlabel("Permeability [D]")
Out[18]:
In [21]:
figure(figsize=(4,3))
semilogy(df_geo[:k_m2]*1e12, df_geo[:COP], "o", alpha = 0.5)
semilogy([0,1.1], [0.8, 0.8], "--k", label = "methane")
semilogy([0,1.1], [1.0, 1.0], "-r", label = "heat pump")
semilogy([0,1.1], [2.0, 2.0], ":g", label = "geo heat pump")
axis([0,1.05,0,40])
legend()
ylabel("COP")
xlabel("Permeability [D]")
savefig("figs/COP_alternatives.png")