-
Notifications
You must be signed in to change notification settings - Fork 0
/
ProductionModel(CCGT-ISCC).R
282 lines (253 loc) · 10 KB
/
ProductionModel(CCGT-ISCC).R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
## load of required libraries
library(raster)
library(solaR)
library(financial)
## introduce the working directory
dir<-as.character('/Users/usuario/Dropbox/Desertec/DocumentosActivos/CCGTAlgeria')
setwd(dir)
## definition of the daily temporal series. In this case 2005 daily time series (dTs)
## is defined with fBTd (solaR function)
dTs<-fBTd(mode='serie',start='01-01-2005',end='31-12-2005',format='%d-%m-%Y')
## definition of the intradaily temporal series (iTs) with as.POSIXct
iTs<-seq(as.POSIXct('2005-01-01 00:30:00'),as.POSIXct('2005-12-31 23:30:00'),by='hour')
## definition of month, day,hour, julian day and julian hour
dom<-c(31,28,31,30,31,30,31,31,30,31,30,31)
month<-rep(seq(1:12),dom*24)
df<-lapply(dom,function(x){t<-as.numeric(seq(1:x))})
day<-rep(do.call(c,df),each=24)
hour<-rep(1:24,365)
julianday<-rep(seq(1:365),each=24)
julianhour<-seq(1:8760)
## Solar field specifications
## reflectivity
LS3reflec<-0.94
## transmisivity
LS3trans<-0.955
## interceptation factor
LS3intf<-0.997
## absortivity
LS3abs<-0.955
## peak efficiency
nopico<-LS3reflec*LS3trans*LS3intf*LS3abs
LS3width<-5.76
LS3focdis<-1.71
LS3length<-99
LS3oparea<-545
mirrorclean<-0.98
##fluid inlet temperature
HTFin<-293
HTFout<-393
HTFm<-343
##other data
PCI<-39.900 ## PCI of natural gas [kJ/kg]
## Supposing cummulative efficiency curves: GTrTaRHP (power)
GTrTaRHP<-function(Ta,RH,P,Pb){
efTa<-(-0.5024*Ta+107.536)/100
efRH<-((1.05/90)*RH+99.3)/100
efP<-((27/25)*((P/Pb)*100-100)+100)/100
ef<-efTa*efRH*efP
ef}
## Supposing cummulative efficiency curves: STrTaRHP (power)
STrTaRHP<-function(Ta){
efTa<-(6e-4*Ta^2-0.1579*Ta+102.24)/100
ef<-efTa
ef}
## definition of stations
ccgt<-as.data.frame(read.csv2('CCGTv1.csv',sep=';',dec=',',header=TRUE),sep=';',header=TRUE,colClasses='numeric',dec='.')
est<-as.character(ccgt$Name[which(ccgt$Land==1)])
###
## definition of lapply function to be applied along 'est'
estccgt<-est[which(ccgt$Cycle=='CCGT')]
for(n_loops in 1:50){
modProd<-lapply(estccgt,function(x){
## introduction of coordinates of the CCGT studied
###
Lat<-as.numeric(as.character(ccgt$Latitude))[which(ccgt$Land==1)][which(est==x)]
Lon<-as.numeric(as.character(ccgt$Longitude))[which(ccgt$Land==1)][which(est==x)]
## introduction of meteo data
setwd('/Users/usuario/Dropbox/Desertec/DocumentosActivos/meteo')
meteo<-as.data.frame(read.csv2(paste0('est',which(est==x),'.csv'),sep=';',header=TRUE,colClasses='numeric',dec='.'))
names(meteo)<-c('y','m','d','h','ghi','TempMed','dni','HumMed','P')
## Apparent movement of the Sun from the Earth: incidence angle
sol<-calcSol(Lat,local2Solar(iTs,Lon),sample='hour',EoT=TRUE,method='michalsky')
incang<-r2d(acos(as.numeric(fTheta(sol,modeTrk='horiz')$cosTheta)[25:8784]))
incang[which(is.na(incang))]<-0
## Definition of modificator of incidence angle as per Gonzalez.Zarza.ea2001
which(incang>80)
## No angle higher than 80?, which would mean modincang<-0
modincang<-1-2.23073e-4*(incang)-1.1e-4*(incang^2)+3.18596e-6*(incang^3)-4.85509e-8*(incang^4)
modincang[which(incang>80)]<-0
## Optical efficiency
optef<-modincang*nopico
## Loss area per collector
Lossarea<-LS3width*(LS3focdis+((LS3focdis*(LS3width^2))/(48*(LS3focdis^2))))*tan(d2r(incang))
## Heat loss from collector to environment (kWth per collector)
Pcolenv<-(0.00154*(HTFm-meteo$TempMed)^2+0.2021*(HTFm-meteo$TempMed)-24.899+
((0.00036*(HTFm-meteo$TempMed)^2+0.2029*(HTFm-meteo$TempMed)+24.899)*(meteo$dni/900)*
cos(d2r(incang))))*LS3length/1000
## Potential thermal Power (kWth per collector)
Psuncol<-LS3oparea*meteo$dni*cos(d2r(incang))/1000
## Thermal power from collector to fluid
Pcolfluid<-(LS3oparea-Lossarea)*meteo$dni*cos(d2r(incang))*optef*mirrorclean/1000-Pcolenv
Pcolfluid[which(Pcolfluid<0)]<-0
Psolar_th<-n_loops*4*Pcolfluid ## 4 Solar Collection Assemblies per loop
## we consider a 34.25% steam turbine efficiency
ef_st<-0.3425
ef_hrsg<-0.95
Psolar_el<-Psolar_th*ef_st*ef_hrsg/1000 ## in MW
## remove NA
Psolar_el[which(is.na(Psolar_el))]<-0
## Annual thermal power from collector to fluid kWhth
Pcolannual<-sum(Pcolfluid,na.rm=TRUE)
## Losses of power of combined cycle, gas turbine and steam turbine
## Definition of inputs
Ta<-meteo$TempMed
RH<-meteo$HumMed
P<-meteo$P
DNI<-meteo$dni
## base pressure
Pb<-((-27/2400)*ccgt$Elevation[which(ccgt$Name==x)]+100)/100*1013
## gas turbine related to air temperature
eGTrTa<-function(x){ef<-(-0.002*x^2-0.1237*x+102.28)}
NG_cons_rate<-(20/260)*as.numeric(as.character(ccgt$TG[which(ccgt$Name==x)]))
efGTrTa<-eGTrTa(Ta)
### SCENARIO 1: solar boosting mode
## a) Natural gas consumption scenario 1
NG_cons_sc1<-NG_cons_rate*efGTrTa/100
## b) GT Power output scenario 1
Pgas_sc1<-GTrTaRHP(Ta,RH,P,Pb)*as.numeric(as.character(ccgt$TG[which(ccgt$Name==x)]))
## c) Integrable solar power scenario 1
Psol_int_sc1<-(-STrTaRHP(Ta)+1)*as.numeric(as.character(ccgt$TV[which(ccgt$Name==x)]))
## remove negative values in Psol_int_sc1
Psol_int_sc1[which(Psol_int_sc1<0)]<-0
## from integrable to energy integrated
Psol_int_sc1i<-lapply(c(1:8760),function(t){
if(Psolar_el[t]>Psol_int_sc1[t]){Psol_int_sc1[t]<-Psol_int_sc1[t]}else{Psol_int_sc1[t]<-Psolar_el[t]}
Total<-Psol_int_sc1[t]})
Psol_int_sc1<-do.call(c,Psol_int_sc1i)
## e) Dumping scenario 1
Dumping_sc1<-Psolar_el-Psol_int_sc1
Dumping_sc1[Dumping_sc1<0]<-0
## d) ST Power output scenario 1
Pst_sc1<-STrTaRHP(Ta)*as.numeric(as.character(ccgt$TV[which(ccgt$Name==x)]))+Psol_int_sc1
## d) Total CCGT power output scenario 1
Pccgt_sc1<-Pst_sc1+Pgas_sc1
## f) Overall efficiency
ef_ccgt_sc1<-Pccgt_sc1/(NG_cons_sc1*PCI)
## Electrical Power
Aux_con<-1.1
ef_sf<-0.52
Rad<-600
Pe_fixed<-n_loops*LS3oparea*4*Rad*ef_st*ef_sf/(Aux_con*1000) # in KW
# conversion price from CCGT to ISCC
p_conv<-3392 # USD/kWe
Cost_fixed<-p_conv*Pe_fixed # cost of hybridation in dollars
## Financial Model
## discount rate
dr<-0.04
infl<-0.02
years<-25
## Positive cash flows
## electricity price for <5% solar share
ep<-0.027214 #USD/kWh
infla<-lapply(c(1:(years-1)),function(x){infl<-(1+infl)^x}); inflation<-do.call(c,infla);inflation<-c(1,inflation)
CF_sell<-sum(Psol_int_sc1*ep*1000)*inflation
## Negative cash flows
CF_fixed<--49.9*Pe_fixed*inflation
CF_var<- -0.7*sum(Psolar_el)*inflation
# Total cash flows
CF0<-Cost_fixed
CF<-CF_sell+CF_fixed+CF_var
residual_cost<-0.1
CF[25]<-CF[25]+residual_cost*CF0
# Residual cost
npv<-sum( do.call(c,lapply(c(1:years),function(t){a<-CF[t]/((1+dr)^t)})))-CF0
CF<-c(-CF0,CF)
irr<-cf(CF,dr)$irr[1]
npv<-rep(npv,8760)
irr<-rep(irr,8760)
# levelized cost of electricity
num<-lapply(c(1:years),function(t){(CF_fixed[t]+CF_var[t])/((1+dr)^t)});num<-do.call(c,num)
num[25]<-CF_fixed[25]+CF_var[25]+residual_cost*CF0
num<--CF0-sum(num)
den<-lapply(c(1:years),function(t){sum(Psol_int_sc1)*1000/((1+dr)^t)});den<-do.call(c,den);den<-sum(den)
lcoe<-num/den
lcoe<-rep(lcoe,8760)
### SCENARIO 0: CCGT
## a) Natural gas consumption
NG_cons_sc0<-NG_cons_rate*efGTrTa/100
## b) GT power output scenario 0
Pgas_sc0<-GTrTaRHP(Ta,RH,P,Pb)*as.numeric(as.character(ccgt$TG[which(ccgt$Name==x)]))
## c) ST power output sceneario 0
Pst_sc0<-STrTaRHP(Ta)*as.numeric(as.character(ccgt$TV[which(ccgt$Name==x)]))
## d) CCGT power output scenario 0
Pccgt_sc0<-Pgas_sc0+Pst_sc0
## e) Solar power integration
Psol_int_sc0<-rep(0,8760)
## g) Overall efficiency
ef_ccgt_sc0<-Pccgt_sc0/(NG_cons_sc0*PCI)
Total<-data.frame(iTs,month,day,julianday,hour,julianhour,DNI,Ta,RH,P,
incang,modincang,optef,Lossarea,Pcolenv,Psuncol,Pcolfluid,Psolar_th,Psolar_el,efGTrTa,
NG_cons_sc1,Pgas_sc1,Psol_int_sc1,Dumping_sc1,Pccgt_sc1,Pst_sc1,ef_ccgt_sc1,
NG_cons_sc0,Pgas_sc0,Pst_sc0,Pccgt_sc0,Psol_int_sc0,ef_ccgt_sc0,npv,irr,lcoe)
Total})
old<-getwd()
setwd('/Users/usuario/Dropbox/Desertec/DocumentosActivos/Resultados/CCGT-ISCC')
save(modProd,file=paste('Todo',n_loops,'.RData',sep=''))
setwd(old)
}
setwd('/Users/usuario/Dropbox/Desertec/DocumentosActivos/Resultados/CCGT-ISCC')
load('Todo1.RData')
horas_operacion<-seq(11,17,1)
horas<-lapply(horas_operacion,function(h){
which(modProd[[1]]$hour==h)})
horas<-sort(do.call(c,horas))
setwd('/Users/usuario/Dropbox/Desertec/DocumentosActivos/Resultados/CCGT-ISCC')
Annual<-lapply(c(1:50),function(x){
load(paste('Todo',x,'.RData',sep=''))
t<- lapply(modProd,function(x){
Pccgt_sc0<-sum(x$Pccgt_sc0[horas])
Pccgt_sc1<-sum(x$Pgas_sc1[horas])+sum(x$Pst_sc1[horas])
Psol_int_sc1<-sum(x$Psol_int_sc1[horas])
Dumping_sc1<-sum(x$Dumping_sc1[horas])
NG_cons_sc0<-sum(x$NG_cons_sc0[horas])
Pgas_sc1<-sum(x$Pgas_sc1[horas])
Pst_sc1<-sum(x$Pst_sc1[horas])
Pst_sc0<-sum(x$Pst_sc0[horas])
Pgt_sc0<-sum(x$Pgas_sc0[horas])
solar_share<-(Psol_int_sc1/Pccgt_sc1)*100
irr<-x$irr[1]
lcoe<-x$lcoe[1]
Total<-c(Pccgt_sc0,Pst_sc0,Pgt_sc0,Pccgt_sc1,Psol_int_sc1,Dumping_sc1,NG_cons_sc0,Pgas_sc1,Pst_sc1,solar_share,irr,lcoe)
})
Total<-data.frame(do.call(rbind,t))
names(Total)<-c('Pccgt_sc0','Pst_sc0','Pgt_sc0','Pccgt_sc1','Psol_int_sc1',
'Dumping_sc1','NG_cons_sc0','Pgas_sc1','Pst_sc1','solar_share','irr','lcoe')
return(Total)
})
save(Annual,file='Annual.RData')
## Selection of the number of loops for which dumping is greater than the solar integration
K<-lapply(c(1:3),function(st_ccgt){
m<-lapply(c(1:50),function(nloop){t<-Annual[[nloop]][st_ccgt,]$Dumping_sc1>Annual[[nloop]][st_ccgt,]$Psol_int_sc1;t})
m<-do.call(c,m)
m<-which(m==TRUE)[1]-1
m<-data.frame(Annual[[m]][st_ccgt,],m)})
## Units
KK<-lapply(K,function(x){
m<-data.frame(x[,c(1,2,3,4,5,6)]/1000,x[,10],x[,12]*100,x[,11],x[,13])
names(m)<-c('Pccgt','Pst','Pgt','Piscc','Psol','Dumping','sol_share','lcoe','irr','nloop')
m
})
est<-c(1,2,3)
R_ccgt<-data.frame(est,do.call(rbind,KK))
print(xtable(R_ccgt,digits=c(1,0,1,1,1,1,1,1,2,2,2,0)),include.rownames=FALSE)
load('Annual.RData')
LCOE<-do.call(rbind,Annual)$lcoe
IRR<-do.call(rbind,Annual)$irr/100
Station<-rep(c('Station 1','Station 2','Station 3'),50)
Loop<-rep(c(1:50),each=3)
LCOE<-data.frame(LCOE,IRR,Station,Loop)
LCOE$Station<-as.factor(LCOE$Station)
trellis.device(pdf, file='LCOEccgt-iscc.pdf')
xyplot(LCOE~Loop,groups=Station,data=LCOE,auto.key=list(space='right', points=TRUE, lines=FALSE))
dev.off()