Carbon Capture

Assigning some numbers:

electricityCost<-0.0477 #$/kWh
span<-20 #yr
componentcosts<-"C:\\Users\\craw038\\Documents\\CarbonCaptureCellCosts.xlsx"%>%
  read_excel()

rownames(componentcosts)<-componentcosts$Component

aspectratio<-1
membraneX<-0.038#m extra length for membrane
frameXW<-0.14 #m extra length for pvc frame
frameXL<-0.05 #m extra length for pvc frame
cdepth<-0.002 #m
cwidth<-0.005 #m


CO2MassRate<-600*24 #kg/day
MMCarbon<-12+16*2 #g/mol CO2
CO2MolarRate<-CO2MassRate*1000/MMCarbon #mol/day
CaptEff<-0.75

Faraday<-26.8 #Ah/mol
CurrentRequired<-CO2MolarRate/24*Faraday/CaptEff #A

CO3Concentration<-2.5/1000 #mol/L
Utilisation<-1


flow<-1/(Faraday*CO3Concentration/1000)/60/Utilisation #ml/min-A
visc<-9e-4 #Pa-s

Voltage<-0.62 #V
j<-(50/1000)/(5/100/100) #A/m^2

Defining Cost function

findCost<-function(Cells,Stacks){
area<-CurrentRequired/Cells/Stacks/j #m^2

q<-CurrentRequired/Cells/Stacks*flow*1e-6/60 #m^3/s

Loss<-function(area,q,Vd){
  
  Dsp<-0.009 #m
  Lsp<-0.1 #m
    Dmp<-1.75*2.54/100 #m
    k<-0.0015/1000 #m roughness
  
  resistivity<-1/5 #S/m
  length<-sqrt(area/aspectratio)
  carea<-Dsp^2*pi/4
  Rs<-resistivity*Lsp/carea
  Rm<-resistivity*Dsp/((Dmp)^2*3.14/4) 
  Ik<-Vd/Rm
  r <- sqrt(Rs/Rm)
  data.frame(cellno=(-(Cells-2)/2):((Cells-2)/2))%>%
    mutate(im=Ik*(1-(cosh(cellno/r)/cosh(0.5*Cells/r))),
           is=(Ik/r)*(sinh(cellno/r)/cosh(0.5*Cells/r)),
           powerLoss=4*Rm*im^2+4*Rs*is^2)%>%
    .$powerLoss%>%mean->shuntLoss
  
  if(visc==0){return(shuntLoss)}
  
  width = sqrt(area*aspectratio)
  height = width/aspectratio
  #Dh=2*cdepth*cwidth/(cdepth+cwidth)
  #carea=cdepth*cwidth
  K = 3.982E-10# m^2
  thickness = 0.001 #m
  #segments <- width/channelsep
  
  rho<-1000 #kg/m^3
  
  ReSP<-q/(Dsp^2*pi/4)*Dsp*rho/visc
  
  froot<-function(f){
    
    if(f<=0){return(NA)}
    1/sqrt(f)+2*log10(k/Dsp/3.7+2.51/ReSP/sqrt(f))
  }
  
  fsp<-uniroot(froot,c(100,1e-10))$root
  
  Lmp<-7 #m
  ReMP<-Cells*q/(Dmp^2*pi/4)*Dmp*rho/visc/2
  
  froot<-function(f){
    
    if(f<=0){return(NA)}
    1/sqrt(f)+2*log10(k/Dmp/3.7+2.51/ReMP/sqrt(f))
  }
  
  fmp<-uniroot(froot,c(100,1e-10))$root
  
  
  DhCell<-2*thickness*width/(thickness+width)
  
  Pdrop<-32*height*q*visc/(thickness*width*DhCell*DhCell)+
    2*fsp*rho*(q/(Dsp^2*pi/4))^2/Dsp/2*Lsp+
    2*fmp*rho*(q*Cells/2/(Dmp^2*pi/4))^2/Dmp/2*Lmp
  
  PEff<-0.67
  powerLoss<-2*Pdrop*q/PEff
  
  return(powerLoss+shuntLoss)
  
}

PowerRequired<-Loss(area,q,Voltage)+area*j*Voltage #W


EnergyCost<-data.frame(Component=c("Electricity","Electricity Loss"),
                       Cost=c(area*j*Voltage,Loss(area,q,Voltage))*Cells*Stacks/1000*electricityCost/(CO2MassRate/24))

length<-sqrt(area*aspectratio)
width<-sqrt(area/aspectratio)

q<-CurrentRequired*flow #ml/min
GPM<-q/3785.41
PumpLoss<-(PowerRequired-area*j*Voltage)*Cells*Stacks/1000/2
#print(PumpLoss)

FeltCost<-componentcosts["Carbon Cloth",]%>%
  mutate(Cost=Cost*area*2*Cells*Stacks)

NfCost<-componentcosts["Membrane",]%>%
  mutate(Cost=2*Cost*(length+membraneX)*(membraneX+width)*Cells*Stacks)


PtCost<-componentcosts["Pt",]%>%
  mutate(Cost=1*Cost*area*Cells*Stacks)

BipolarCost<-componentcosts["Carbon Plate",]%>%
  mutate(Cost=Cost*area*Cells*Stacks)

PVCCost<-componentcosts["PVC Frame",]%>%
  mutate(Cost=Cost*(length+frameXW)*(width+frameXL)*2*Cells*Stacks)

PumpCost<-componentcosts["Pump Cost per GPM",]%>%
  mutate(Component="Pumps",
         Cost=(Cost*GPM*2+componentcosts["Pump Base Cost",]$Cost*2))

PumpCost<-componentcosts["Pump Cost per kW",]%>%
  mutate(Component="Pumps",
         Cost=(Cost*PumpLoss*2+componentcosts["Pump Base Cost2",]$Cost*2))

PCSCost<-componentcosts["PCS",]%>%
  mutate(Cost=Cost*PowerRequired*Cells*Stacks)

AlCosts<-componentcosts["Al Endplate",]%>%
  mutate(Cost=Cost*2*Stacks)

bind_rows(EnergyCost,
          PtCost,
          FeltCost,
          NfCost,
          BipolarCost,
          PVCCost,
          PumpCost,
          PCSCost,
          AlCosts)%>%
  return
}

Levelised Cost Multipliers:

"C:\\Users\\craw038\\Documents\\LevelisedCost.csv"%>%
  read.csv(row.names = NULL)->X

X%>%lm(data=.,`ï..Levelised.Cost`~.+0)%>%summary%>%coef%>%.[,"Estimate"]->LCOEMult

LCOEMultiplier<-bind_rows(data.frame(Component=c("Electricity",
                             "Electricity Loss"),
                           Mult=CO2MassRate*365.25*LCOEMult["Electricity"]),
                          data.frame(Component=c("PVC Frame",
                                                 "Pumps",
                                                 "PCS",
                                                 "AI Endplate"),
                                     Mult=LCOEMult["Capital"]),
                          data.frame(Component=c("Membrane",
                                                 "Pt",
                                                 "Carbon Cloth",
                                                 "Carbon Plate"),
                                     Mult=LCOEMult["CapitalReplace"]))

Check various electrochemical performance scenarios:

ElectroChemScenarios<-data.frame(V=c(0.62,0.88,0.62,0.93),
                                 j=c(5,10,10,50)*100*100/1000)

ElectroChemScenarios<-bind_rows(ElectroChemScenarios%>%mutate(Pumped=TRUE),
                                ElectroChemScenarios%>%mutate(Pumped=FALSE))

ElectroChemResults<-NULL
for(i in 1:nrow(ElectroChemScenarios)){
  
  Voltage<-ElectroChemScenarios$V[i] #V
  j<-ElectroChemScenarios$j[i] #A/m^2
  visc <- 0.9E-03*ElectroChemScenarios$Pumped[i] # p
  
  CO2MassRate<-24*ifelse(ElectroChemScenarios$Pumped[i],750,450) #kg/day
  MMCarbon<-12+16*2 #g/mol CO2
  CO2MolarRate<-CO2MassRate*1000/MMCarbon #mol/day
  CaptEff<-0.75
  
  Faraday<-26.8 #Ah/mol
  CurrentRequired<-CO2MolarRate/24*Faraday/CaptEff #A
  
  CO3Concentration<-2.5/1000 #mol/L
  Utilisation<-1
  electricityCost<-0.0477 #$/kWh
  
  flow<-1/(Faraday*CO3Concentration/1000)/60/Utilisation #ml/min-A
  
  
  findCost(15,800)%>%
    mutate(V=Voltage,
           j=j,
           Pumped=ElectroChemScenarios$Pumped[i],
           Scenario=paste0(Voltage," V\n", j/10, " mA/cm^2"),
           CO2perHour=CO2MassRate/24)%>%
    bind_rows(ElectroChemResults)->ElectroChemResults
  
}

ElectroChemResults%>%
  filter(!(Component=="Pumps" & !Pumped))%>%
  merge(LCOEMultiplier)%>%
  mutate(LevlisedCost=Cost*Mult*ifelse(!grepl("Electricity",Component),600/CO2perHour,1),
         Scenario2=ifelse(Pumped,"Pumped","Effulent"),
         Component=factor(Component, levels(factor(Component))[rev(c(5,1,2,6:9,3,4))]))%>%
    select(-Cost)%>%
  arrange(Pumped,Scenario,Component)%T>%write.csv("CarbonCaptureCostBase.csv")%>%
  ggplot(aes(x=Scenario,y=LevlisedCost,fill=Component))+
  geom_bar(stat="identity", colour="black",size=1)+theme_bw()+ylab("Levelized Cost $/kg")+facet_wrap(.~Scenario2)+scale_fill_brewer(palette = "Set1")+
    theme(strip.background = element_blank(),
          axis.title = element_text(face="bold",size=16),
          axis.text = element_text(face="bold",size=10,colour="black"),
          strip.text =  element_text(face="bold",size=16,colour="black"),
          legend.text =   element_text(face="bold",size=14,colour="black"),
          legend.title =  element_text(face="bold",size=14,colour="black"))+
    scale_y_continuous(expand = expand_scale(mult = c(0,0.05), add = 0), limits = c(0, NA))

ElectroChemResults%>%
  filter(!(Component=="Pumps" & !Pumped))%>%
  merge(LCOEMultiplier)%>%
  mutate(LevlisedCost=Cost*Mult*ifelse(!grepl("Electricity",Component),600/CO2perHour,1),
         Scenario2=ifelse(Pumped,"Pumped","Effulent"),
         Component=factor(Component, levels(factor(Component))[rev(c(5,1,2,6:9,3,4))]))%>%
  filter(!grepl("Electricity",Component))%>%
    select(-LevlisedCost)%>%
  arrange(Pumped,Scenario,Cost)%T>%write.csv("CarbonCaptureCostBaseCapital.csv")%>%
  ggplot(aes(x=Scenario,y=Cost,fill=Component))+
  geom_bar(stat="identity", colour="black",size=1)+theme_bw()+ylab("Levelized Cost $/kg")+facet_wrap(.~Scenario2)+scale_fill_brewer(palette = "Set1")+
    theme(strip.background = element_blank(),
          axis.title = element_text(face="bold",size=16),
          axis.text = element_text(face="bold",size=10,colour="black"),
          strip.text =  element_text(face="bold",size=16,colour="black"),
          legend.text =   element_text(face="bold",size=14,colour="black"),
          legend.title =  element_text(face="bold",size=14,colour="black"))+
    scale_y_continuous(expand = expand_scale(mult = c(0,0.05), add = 0), limits = c(0, NA))

Load in electricity price data:

MarketData<-read_feather("C:\\Users\\craw038\\Documents\\RedoxCellCost\\RegionalElectricityPrice.feather")

Scenarios2<-NULL
for(i in 1:24){
  
  MarketData%>%
    arrange(Region,Day,Cost)%>%
    group_by(Region,Day)%>%
    mutate(Rank=seq_along(Cost),
           Count=n(),
           Cost=Cost/1000)%>%
    filter(Count==24,
           Rank<=i)%>%
    group_by(Region)%>%
    summarise(MeanCost=mean(Cost))%>%
    mutate(Hours=i)%>%
    bind_rows(Scenarios2)->Scenarios2
}

Scenarios2<-bind_rows(Scenarios2%>%mutate(Pumped=TRUE),Scenarios2%>%mutate(Pumped=FALSE))

Scenarios2%>%ggplot(aes(x=Hours,y=MeanCost,colour=Region))+geom_line(size=1.2)+
  theme_bw()+
  xlab("Hours per Day Operated")+
  ylab("Electricity Cost ($/kWh)")+
  scale_color_brewer(palette="Set1")+
    theme(strip.background = element_blank(),
          axis.title = element_text(face="bold",size=16),
          axis.text = element_text(face="bold",size=10,colour="black"),
          strip.text =  element_text(face="bold",size=16,colour="black"),
          legend.text =   element_text(face="bold",size=14,colour="black"),
          legend.title =  element_text(face="bold",size=14,colour="black"))+
    scale_y_continuous(expand = expand_scale(mult = c(0,0.05), add = 0), limits = c(0, NA))

Check various electricity price scenarios:

Costs2<-NULL
for(i in 1:nrow(Scenarios2)){
  
  
Voltage<-0.62 #V
j<-(50/1000)/(5/100/100) #A/m^2
  
  CO2MassRate<-ifelse(Scenarios2$Pumped[i],750,450)*24*(24/Scenarios2$Hours[i]) #kg/day
  MMCarbon<-12+16*2 #g/mol CO2
  CO2MolarRate<-CO2MassRate*1000/MMCarbon #mol/day
  CaptEff<-0.75
  
  Faraday<-26.8 #Ah/mol
  CurrentRequired<-CO2MolarRate/24*Faraday/CaptEff #A
  
  CO3Concentration<-2.5/1000 #mol/L
  Utilisation<-1

  flow<-1/(Faraday*CO3Concentration/1000)/60/Utilisation #ml/min-A

  electricityCost<-Scenarios2$MeanCost[i]
  
  
Costs<-findCost(15,800)

Costs2<-Costs%>%mutate(Region=Scenarios2$Region[i],
                       Hours=Scenarios2$Hours[i],
                       Pumped=Scenarios2$Pumped[i],
                        CO2perHour=ifelse(Scenarios2$Pumped[i],750,450))%>%bind_rows(Costs2)

}

Costs2%>%
  filter(!(Component=="Pumps" & !Pumped))%>%
  merge(LCOEMultiplier)%>%
  mutate(LevlisedCost=Cost*Mult*ifelse(!grepl("Electricity",Component),600/CO2perHour,1),
         Scenario2=ifelse(Pumped,"Pumped","Effulent"),
         Component=factor(Component, levels(factor(Component))[rev(c(5,1,2,6:9,3,4))]))%>%
  select(-Cost)%>%
  arrange(Region,Scenario2,Hours,Component)%T>%write.csv("CarbonCaptureCostsbyRegion.csv")%>%
  group_by(Region,Scenario2,Hours)%>%
  mutate(TotalCost=sum(LevlisedCost))%>%
  group_by(Region,Scenario2)%>%#filter(Scenario2=="Pumped",Region=="PNW")%>%
  filter(Hours==24 | TotalCost==min(TotalCost))%>%
  mutate(Scenario=ifelse(Hours==24 & !((Hours%>%unique%>%sum)==24),"Non-Opt","Opt"))%>%
  #mutate(Cost=Cost/ifelse(grepl("Electricity",Component),1,CO2MassRate*365.25*20))%>%
  ggplot(aes(x=Scenario,y=LevlisedCost,fill=Component))+
  geom_bar(stat="identity", colour="black",size=1)+theme_bw()+ylab("Levelized Cost $/kg")+facet_grid(Scenario2~Region)+scale_fill_brewer(palette = "Set1")+
    theme(strip.background = element_blank(),
          axis.title = element_text(face="bold",size=16),
          axis.text = element_text(face="bold",size=10,colour="black"),
          strip.text =  element_text(face="bold",size=16,colour="black"),
          legend.text =   element_text(face="bold",size=14,colour="black"),
          legend.title =  element_text(face="bold",size=14,colour="black"))+
    scale_y_continuous(expand = expand_scale(mult = c(0,0.05), add = 0), limits = c(0, NA))

Using NF212 and 10x Pt cost:

ElectroChemScenarios<-data.frame(V=c(0.62,0.88,0.62,0.93),
                                 j=c(5,10,10,50)*100*100/1000)

ElectroChemScenarios<-bind_rows(ElectroChemScenarios%>%mutate(Pumped=TRUE),
                                ElectroChemScenarios%>%mutate(Pumped=FALSE))

ElectroChemResults<-NULL
for(i in 1:nrow(ElectroChemScenarios)){
  
  Voltage<-ElectroChemScenarios$V[i] #V
  j<-ElectroChemScenarios$j[i] #A/m^2
  visc <- 0.9E-03*ElectroChemScenarios$Pumped[i] # p
  
  CO2MassRate<-24*ifelse(ElectroChemScenarios$Pumped[i],750,450) #kg/day
  MMCarbon<-12+16*2 #g/mol CO2
  CO2MolarRate<-CO2MassRate*1000/MMCarbon #mol/day
  CaptEff<-0.75
  
  Faraday<-26.8 #Ah/mol
  CurrentRequired<-CO2MolarRate/24*Faraday/CaptEff #A
  
  CO3Concentration<-2.5/1000 #mol/L
  Utilisation<-1
  electricityCost<-0.0477 #$/kWh
  
  flow<-1/(Faraday*CO3Concentration/1000)/60/Utilisation #ml/min-A
  
  
  findCost(15,800)%>%
    mutate(V=Voltage,
           j=j,
           Pumped=ElectroChemScenarios$Pumped[i],
           Scenario=paste0(Voltage," V\n", j/10, " mA/cm²"),
           CO2perHour=CO2MassRate/24)%>%
    bind_rows(ElectroChemResults)->ElectroChemResults
  
}

(ElectroChemResults%>%
  filter(!(Component=="Pumps" & !Pumped))%>%
  merge(LCOEMultiplier)%>%
  mutate(Cost=Cost*ifelse(Component=="Membrane",95/32,1),
         Cost=Cost*ifelse(Component=="Pt",10,1))%>%
  mutate(LevlisedCost=Cost*Mult*ifelse(!grepl("Electricity",Component),600/CO2perHour,1),
         Scenario2=ifelse(Pumped,"Pumped","Effulent"),
         Component=factor(Component, levels(factor(Component))[rev(c(5,1,2,6:9,3,4))]))%>%
    select(-Cost)%>%
  arrange(Pumped,Scenario,Component)%T>%write.csv("CarbonCaptureCostNf212Pt10x.csv")%>%
  ggplot(aes(x=Scenario,y=LevlisedCost,fill=Component))+
  geom_bar(stat="identity", colour="black",size=1)+theme_bw()+ylab("Levelized Cost $/kg")+facet_wrap(.~Scenario2)+scale_fill_brewer(palette = "Set1")+
    theme(strip.background = element_blank(),
          axis.title = element_text(face="bold",size=16),
          axis.text = element_text(face="bold",size=10,colour="black"),
          strip.text =  element_text(face="bold",size=16,colour="black"),
          legend.text =   element_text(face="bold",size=14,colour="black"),
          legend.title =  element_text(face="bold",size=14,colour="black"))+
    scale_y_continuous(expand = expand_scale(mult = c(0,0.05), add = 0), limits = c(0, NA)))%>%print

ElectroChemResults%>%
  filter(!(Component=="Pumps" & !Pumped))%>%
  merge(LCOEMultiplier)%>%
  mutate(LevlisedCost=Cost*Mult*ifelse(!grepl("Electricity",Component),600/CO2perHour,1),
         Scenario2=ifelse(Pumped,"Pumped","Effulent"),
         Component=factor(Component, levels(factor(Component))[rev(c(5,1,2,6:9,3,4))]))%>%
   mutate(Cost=Cost*ifelse(Component=="Membrane",95/32,1),
         Cost=Cost*ifelse(Component=="Pt",10,1))%>%
  filter(!grepl("Electricity",Component))%>%
    select(-LevlisedCost)%>%
  arrange(Pumped,Scenario,Cost)%T>%write.csv("CarbonCaptureCostNf212Pt10xCapital.csv")%>%
  ggplot(aes(x=Scenario,y=Cost,fill=Component))+
  geom_bar(stat="identity", colour="black",size=1)+theme_bw()+ylab("Capital Cost $")+facet_wrap(.~Scenario2)+scale_fill_brewer(palette = "Set1")+
    theme(strip.background = element_blank(),
          axis.title = element_text(face="bold",size=16),
          axis.text = element_text(face="bold",size=10,colour="black"),
          strip.text =  element_text(face="bold",size=16,colour="black"),
          legend.text =   element_text(face="bold",size=14,colour="black"),
          legend.title =  element_text(face="bold",size=14,colour="black"))+
    scale_y_continuous(expand = expand_scale(mult = c(0,0.05), add = 0), limits = c(0, NA))

Costs2<-NULL
for(i in 1:nrow(Scenarios2)){
  
  
Voltage<-0.62 #V
j<-(50/1000)/(5/100/100) #A/m^2
  
  CO2MassRate<-ifelse(Scenarios2$Pumped[i],750,450)*24*(24/Scenarios2$Hours[i]) #kg/day
  MMCarbon<-12+16*2 #g/mol CO2
  CO2MolarRate<-CO2MassRate*1000/MMCarbon #mol/day
  CaptEff<-0.75
  
  Faraday<-26.8 #Ah/mol
  CurrentRequired<-CO2MolarRate/24*Faraday/CaptEff #A
  
  CO3Concentration<-2.5/1000 #mol/L
  Utilisation<-1

  flow<-1/(Faraday*CO3Concentration/1000)/60/Utilisation #ml/min-A

  electricityCost<-Scenarios2$MeanCost[i]
  
  
Costs<-findCost(15,800)

Costs2<-Costs%>%mutate(Region=Scenarios2$Region[i],
                       Hours=Scenarios2$Hours[i],
                       Pumped=Scenarios2$Pumped[i],
                        CO2perHour=ifelse(Scenarios2$Pumped[i],750,450))%>%bind_rows(Costs2)

}

Costs2%>%
  filter(!(Component=="Pumps" & !Pumped))%>%
  merge(LCOEMultiplier)%>%
    mutate(Cost=Cost*ifelse(Component=="Membrane",95/32,1),
         Cost=Cost*ifelse(Component=="Pt",10,1))%>%
  mutate(LevlisedCost=Cost*Mult*ifelse(!grepl("Electricity",Component),600/CO2perHour,1),
         Scenario2=ifelse(Pumped,"Pumped","Effulent"),
         Component=factor(Component, levels(factor(Component))[rev(c(5,1,2,6:9,3,4))]))%>%
  select(-Cost)%>%
  arrange(Region,Scenario2,Hours,Component)%T>%write.csv("CarbonCaptureCostsbyRegionNf212Pt10x.csv")%>%
  group_by(Region,Scenario2,Hours)%>%
  mutate(TotalCost=sum(LevlisedCost))%>%
  group_by(Region,Scenario2)%>%#filter(Scenario2=="Pumped",Region=="PNW")%>%
  filter(Hours==24 | TotalCost==min(TotalCost))%>%
  mutate(Scenario=ifelse(Hours==24 & !((Hours%>%unique%>%sum)==24),"Non-Opt","Opt"))%>%
  #mutate(Cost=Cost/ifelse(grepl("Electricity",Component),1,CO2MassRate*365.25*20))%>%
  ggplot(aes(x=Scenario,y=LevlisedCost,fill=Component))+
  geom_bar(stat="identity", colour="black",size=1)+theme_bw()+ylab("Levelized Cost $/kg")+facet_grid(Scenario2~Region)+scale_fill_brewer(palette = "Set1")+
    theme(strip.background = element_blank(),
          axis.title = element_text(face="bold",size=16),
          axis.text = element_text(face="bold",size=10,colour="black"),
          strip.text =  element_text(face="bold",size=16,colour="black"),
          legend.text =   element_text(face="bold",size=14,colour="black"),
          legend.title =  element_text(face="bold",size=14,colour="black"))+
    scale_y_continuous(expand = expand_scale(mult = c(0,0.05), add = 0), limits = c(0, NA))