#!/usr/bin/env python # coding: utf-8 # # Padulles-Amphlett Dynamic Model # ### Version 1.0 # # ## Overview #

# The Padulles dynamic model can predict the transient response of cell voltage, temperature of the cell, hydrogen/oxygen out flow rates and cathode and anode channel temperatures/pressures under sudden change in load current. Hence, a dynamic fuel cell simulation is developed in this model, which incorporates the dynamics of flow and pressure in the anode and cathode channels and mass/ heat transfer transient features in the fuel cell body. #
This model is based on several assumptions: #

    #
  1. The stack is fed with hydrogen and air
  2. #
  3. Cell temperature is stable at all times
  4. #
  5. The ratio of pressures between the interior and exterior of the electrode channels is large
  6. #
  7. The channels that transport gases along the electrodes have a fixed volume
  8. #
  9. Only source of voltage loss is ohmic polarization
  10. #
  11. Nernst equation can be applied too
  12. #
#

#

# This model is an integration of Padulles-Hauer dynamic model with Amphlett static model. The advantage of this dynamic model is using Amphlett equation for simulating the polarization values. Amphlett model as the most complicated and preferable static model, but the most precise. Based on this model, the obtained polarization voltage is identical to the experimental results. #

#
# # #

Fig1. Padulles-Amphlett Dynamic Model Block Diagram

# #
# ## Nernst Voltage # $$E_{Nernst}=N_0\times [E_0+\frac{RT}{2F}ln(\frac{P_{H_2}\times \sqrt{P_{O_2}}}{P_{H_2O}})]$$ # $$P_{H_2}=\frac{\frac{1}{K_{H_2}}}{1+\tau_{H_2}^{(s)}}[(q_{H_2}^{(inlet)}-(2\times K_r \times i)]$$ # $$\frac{q_{H_2}^{(inlet)}}{q_{methanol}}=\frac{CV}{\tau_{1}^{(s)}+(\tau_{2}^{(s)})^2+(\tau_{1}+\tau_{2})^{(s)}+1}$$ # $$P_{O_2}=\frac{\frac{1}{K_{O_2}}}{1+\tau_{O_2}^{(s)}}[(q_{O_2}^{(inlet)}-(K_r \times i)]$$ # $$P_{H_2O}=\frac{\frac{1}{K_{H_2O}}}{1+\tau_{H_2O}^{(s)}}[(q_{H_2O}^{(inlet)}-(2\times K_r \times i)]$$ # $$K_r=\frac{N_0}{4F}$$ # $$q_{O_2}^{(inlet)}=\frac{q_{H_2}^{(inlet)}}{r_{h-o}}$$ # $$q_{H_2O}^{(inlet)}=q_{H_2}^{(inlet)}$$ # In[1]: from opem.Dynamic.Padulles_Amphlett import qH2_Calc,Kr_Calc,qO2_Calc,PH2_Calc,PO2_Calc,PH2O_Calc,Enernst_Calc # In[2]: qH2=qH2_Calc(qMethanol=0.0002,CV=2,t1=2,t2=2) qH2 # In[3]: Kr=Kr_Calc(N0=5) Kr # In[4]: qO2=qO2_Calc(qH2=qH2,rho=1.168) qO2 # In[5]: PH2=PH2_Calc(KH2=0.0000422,tH2=3.37,Kr=Kr,I=1,qH2=qH2) PH2 # In[6]: PO2=PO2_Calc(KO2=0.0000211,tO2=6.74,Kr=Kr,I=1,qO2=qO2) PO2 # In[7]: PH2O=PH2O_Calc(KH2O=0.000007716,tH2O=18.418,Kr=Kr,I=1,qH2O=qH2) PH2O # In[8]: Enernst=Enernst_Calc(E0=0.6,N0=5,T=343, PH2=PH2, PO2=PO2,PH2O=PH2O) Enernst # ## PEM Losses Model # ### Activation # $$\eta_{activation}=\xi_{1}+\xi_{2}T+\xi_{3}T[ln(C_{O_{2}})]+\xi_{4}T[ln(i)]$$ # $$\xi_{1}=-0.948$$ # $$\xi_{2}=0.00286+0.0002\times ln(A)+(4.3\times10^{-5})[ln(C_{H_{2}})]$$ # $$\xi_{3}=7.6\times10^{-5}$$ # $$\xi_{4}=-1.93\times10^{-4}$$ # $$C_{H_{2}}=\frac{P_{H_2}}{1.09\times10^{6}\times exp(\frac{77}{T})}$$ # $$C_{O_{2}}=\frac{P_{O_2}}{5.08\times10^{6}\times exp(\frac{-498}{T})}$$ # In[9]: from opem.Dynamic.Padulles_Amphlett import Eta_Act_Calc Eta_Act=Eta_Act_Calc(T=343,PO2=PO2 , PH2=PH2, i=2, A=50.6) Eta_Act # ### Ohmic # $$\eta_{ohmic}=i(R_{electronic}+R_{Proton})$$ # $$R_{Proton}=\frac{\rho_m\times I}{A}$$ # $$\rho_m=\frac{181.6[1+0.03(\frac{i}{A})+0.062(\frac{T}{303})^2(\frac{i}{A})^{2.5}]}{[\lambda-0.634-3(\frac{i}{A})]exp[4.18(\frac{T-303}{T})]}$$ # * R-Electronic Should be approximately constant over the relatively narrow # temperature range of PEM fuel cell operation. # Therefore, the parameter R-Electronic can be taken as a constant, but # is generally difficult to predict and, therefore, is initially an unknown. # In[10]: from opem.Dynamic.Padulles_Amphlett import Eta_Ohmic_Calc Eta_Ohmic=Eta_Ohmic_Calc(i=2, l=0.0178, A=50.6, T=343, lambda_param=23, R_elec=0) Eta_Ohmic # ### Concentration # $$\eta_{Concentration}=-B\times ln(1-\frac{J}{J_{Max}})$$ # $$J=\frac{i}{A}$$ # $$J_{Max}=\frac{i_L}{A}$$ # In[11]: from opem.Dynamic.Padulles_Amphlett import Eta_Conc_Calc Eta_Conc=Eta_Conc_Calc(i=2, A=50.6, B=0.016, JMax=1.5) Eta_Conc # ### FC Voltage # $$Loss=\eta_{Activation}+\eta_{Ohmic}+\eta_{Concentration}$$ # $$V_{Fuelcell}=E_{Nernst}-N_0\times Loss$$ # In[12]: from opem.Dynamic.Padulles_Amphlett import Loss_Calc,Vcell_Calc Loss=Loss_Calc(Eta_Act=Eta_Act,Eta_Conc=Eta_Conc,Eta_Ohmic=Eta_Ohmic) Loss # In[13]: FC_Voltage=Vcell_Calc(Enernst=Enernst, Loss=Loss, N=5) FC_Voltage # ## Power of PEMFC # $$P=V_{cell}\times i$$ # $$P_{Thermal}=i\times(N \times E_{th}-V_{Stack})$$ # $$E_{th}=\frac{-\Delta H}{nF}=1.23V$$ # In[14]: from opem.Dynamic.Padulles_Amphlett import Power_Calc,Power_Thermal_Calc Power=Power_Calc(Vcell=FC_Voltage,i=2) Power # In[15]: Power_Thermal_Calc(VStack=FC_Voltage,N=5,i=2) # ## Efficiency of PEMFC # $$\eta=\frac{\mu_F\times V_{Fuelcell}}{N_0\times HHV}$$ # In[16]: from opem.Dynamic.Padulles_Amphlett import Efficiency_Calc Efficiency_Calc(Vcell=FC_Voltage,N=5) # ## Linear Approximation # Sometimes quick calculations regarding fuel cell efficiency–power-size relationships need to be made. Linear approximation is a good method to find a rough estimate of the value of polarization function at a particular point. A linear polarization curve has the following form: # $$V_{cell}=V_0-kI$$ # where V0 is the intercept (actual open circuit voltage is always higher) and k is the slope of the curve. # # * Notice : Simple linear regression used for this approximation # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
$$Parameter$$$$Description$$$$Unit$$
$$V_0$$Intercept of the curve obtained by linear approximation$$V$$
$$k$$Slope of the curve obtained by linear approximation$$A^{-1}$$
$$P_{max}$$Maximum power obtained by linear approximation$$W$$
$$V_{FC}|P_{max}$$Cell voltage at maximum power obtained by linear approximation$$V$$
# * Notice : These parameters are only available in HTML report # ## Overall Parameters # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
$$Parameter$$$$Description$$$$Unit$$
$$\eta|P_{Max}$$Cell efficiency at maximum power$$--$$
$$P_{Max}$$Maximum power $$W$$
$$P_{Elec} $$Total electrical power $$W$$
$$P_{Thermal} $$Total thermal power $$W$$
$$V_{FC}|P_{Max}$$Cell voltage at maximum power $$V$$
# * Notice : P(Thermal) & P(Elec) calculated by Simpson's Rule # # * Notice : These parameters are only available in HTML report # ## Full Run # * Run from `i`=0.1 to `i`=75 with `step`=0.1 # In[17]: Test_Vector = { "A": 50.6, "l": 0.0178, "lambda": 23, "JMax": 1.5, "T": 343, "N0": 5, "KO2": 0.0000211, "KH2": 0.0000422, "KH2O": 0.000007716, "tH2": 3.37, "tO2": 6.74, "t1": 2, "t2": 2, "tH2O": 18.418, "B": 0.016, "rho": 1.168, "qMethanol": 0.0002, "CV": 2, "i-start": 0.1, "i-stop": 75, "i-step": 0.1, "Name": "Padulles_Amphlett_Test"} # In[18]: from opem.Dynamic.Padulles_Amphlett import Dynamic_Analysis data=Dynamic_Analysis(InputMethod=Test_Vector,TestMode=True,PrintMode=True,ReportMode=True) # * Notice : "Status", "V0", "K" and "EFF" , new in version 0.8 # In[19]: data_2=Dynamic_Analysis(InputMethod=Test_Vector,TestMode=True,PrintMode=True,ReportMode=False) # * Notice : "PrintMode" & "ReportMode" , new in version 0.5 # In[20]: Dynamic_Analysis(InputMethod={},TestMode=True,PrintMode=False,ReportMode=True) # * Notice: # 1. `TestMode` : Active test mode and get/return data as `dict`, (Default : `False`) # 2. `ReportMode` : Generate reports(`.csv`,`.opem`,`.html`) and print result in console, (Default : `True`) # 3. `PrintMode` : Control printing in console, (Default : `True`) # ## Plot # In[21]: import sys get_ipython().system('{sys.executable} -m pip install matplotlib;') import matplotlib.pyplot as plt # In[22]: def plot_func(x,y,x_label,y_label,color='green',legend=[],multi=False): plt.figure() plt.grid() if multi==True: for index,y_item in enumerate(y): plt.plot(x,y_item,color=color[index]) else: plt.plot(x,y,color=color) plt.xlabel(x_label) plt.ylabel(y_label) if len(legend)!=0: plt.legend(legend) plt.show() # In[23]: plot_func(data["I"],data["P"],"I(A)","Power(W)") # In[24]: plot_func(data["I"],[data["V"],data["VE"]],"I(A)","V(V)",["red","blue"],legend=["Voltage-Stack","Linear-Apx"],multi=True) # In[25]: plot_func(data["I"],data["Ph"],"I(A)","Power-Thermal(W)","blue") # In[26]: plot_func(data["I"],data["EFF"],"I(A)","Efficiency","gray") # In[27]: plot_func(data["I"],data["PO2"],"I(A)","PO2(atm)","orange") # In[28]: plot_func(data["I"],data["PH2"],"I(A)","PH2(atm)","violet") # In[29]: plot_func(data["I"],data["PH2O"],"I(A)","PH2O(atm)","pink") # In[30]: plot_func(data["I"],[data["Eta_Active"],data["Eta_Conc"],data["Eta_Ohmic"]],"I(A)","V(V)",["red","blue","orange"], legend=["Eta_Active","Eta_Conc","Eta_Ohmic"],multi=True) # HTML File # OPEM File # CSV File # ## Parameters # Inputs, Constants & Middle Values # 1. User : User input # 2. System : Simulator calculation (middle value) # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
$$Parameter$$$$Description$$$$Unit$$$$Value$$
$$T$$Fuel cell temperature$$K$$$$User$$
$$N_0$$Number of cells$$--$$$$User$$
$$E_0$$No load voltage$$V$$$$User$$
$$K_{H_2}$$Hydrogen valve constant$$kmol.s^{-1}.atm^{-1}$$$$User$$
$$K_{H_2O}$$Water valve constant$$kmol.s^{-1}.atm^{-1}$$$$User$$
$$K_{O_2}$$Oxygen valve constant$$kmol.s^{-1}.atm^{-1}$$$$User$$
$$\tau_{H_2}^{(s)}$$Hydrogen time constant$$s$$$$User$$
$$\tau_{H_2O}^{(s)}$$Water time constant$$s$$$$User$$
$$\tau_{O_2}^{(s)}$$Oxygen time constant$$s$$$$User$$
$$l$$Membrane thickness$$cm$$$$User$$
$$A$$Active area$$cm^2$$$$User$$
$$\tau_{1}^{(s)}$$Reformer time constant$$s$$$$User$$
$$\tau_{2}^{(s)}$$Reformer time constant$$s$$$$User$$
$$CV$$Conversion factor$$--$$$$User$$
$$B$$An empirical constant # depending on the cell and its # operation state$$V$$$$User$$
$$R_{electronic}$$R-Electronic$$\Omega$$$$User$$
$$\lambda$$An adjustable parameter with a possible minimum value of 14 and a maximum value of 23$$--$$$$User$$
$$J_{Max}$$Maximum current density of the cell$$Acm^{-2}$$$$User$$
$$r_{h-o}$$Hydrogen-Oxygen flow ratio$$--$$$$User$$
$$q_{methanol}$$Molar flow of methanol$$kmol.s^{-1}$$$$User$$
$$i_{start}$$Cell operating current start point$$A$$$$User$$
$$i_{step}$$Cell operating current step$$A$$$$User$$
$$i_{stop}$$Cell operating current end point$$A$$$$User$$
$$P_{H_2}$$Hydrogen partial pressure$$atm$$$$System$$
$$P_{H_2O}$$Water partial pressure$$atm$$$$System$$
$$P_{O_2}$$Oxygen partial pressure$$atm$$$$System$$
$$K_r$$Modeling constant$$kmol.s^{-1}.A^{-1}$$$$System$$
$$q_{O_2}^{(inlet)}$$Molar flow of oxygen$$kmol.s^{-1}$$$$System$$
$$q_{H_2O}^{(inlet)}$$Molar flow of water$$kmol.s^{-1}$$$$System$$
$$q_{H_2}^{(inlet)}$$Molar flow of hydrogen$$kmol.s^{-1}$$$$System$$
$$J$$Actual current density of the cell $$Acm^{-2}$$$$System$$
$$C_{O_2}$$Concentration of oxygen in the catalytic interface of the cathode$$molcm^{-3}$$$$System$$
$$C_{H_2}$$Concentration of hydrogen in the catalytic interface of the anode$$molcm^{-3}$$$$System$$
$$R_{Proton}$$Resistance to proton flow$$\Omega$$$$System$$
$$\xi_2$$Parametric coefficients for cell model$$--$$$$System$$
$$\xi_1$$Parametric coefficients for cell model$$--$$$$-0.948$$
$$\xi_3$$Parametric coefficients for cell model$$--$$$$7.6\times10^{-5}$$
$$\xi_4$$Parametric coefficients for cell model$$--$$$$-1.93\times10^{-4}$$
$$\mu_F$$The fuel utilization$$--$$$$0.95$$
$$HHV$$Higher heating value potential$$V$$$$1.482$$
$$R$$Universal gas constant$$J.kmol^{-1}.K^{-1}$$$$8314.47$$
$$F$$Faraday’s constant$$C.kmol^{-1}$$$$96484600$$
$$E_{th}$$Theoretical potential$$V$$$$1.23$$
# * Notice : $$q_{H_2}=q_{H_2O}$$ # ## Reference #
# 1- J. Padulles, G.W. Ault, J.R. McDonald. 2000. "An integrated SOFC plant dynamic model for power systems # simulation." Journal of Power Sources (Elsevier) 86 (1-2): 495-500. doi:10.1016/S0378-7753(99)00430-9 #
#
# 2- Hauer, K.-H. 2001. "Analysis tool for fuel cell vehicle hardware and software (controls) with an application # to fuel economy comparisons of alternative system designs." Ph.D. dissertation, Transportation Technology # and Policy, University of California Davis. #
#
# 3- J. C. Amphlett, R. M. Baumert, R. F. Mann, B. A. Peppley, and P. R. Roberge. 1995. "Performance Modeling # of the Ballard Mark IV Solid Polymer Electrolyte Fuel Cell." J. Electrochem. Soc. (The Electrochemical Society, # Inc.) 142 (1): 9-15. doi: 10.1149/1.2043959. #