Merge pull request #4 from c-randall/surf_chem

update for paper submission

Looping_Runners/w_Pt_runner.py

@@ -26,24 +26,24 @@
 plt.rcParams.update({'font.size': font_sz})
 
 " Loop through Pt loading polarization curves: "
-#folder_nm_generic = 'cs_Pt_loop_' # pre-fix to all folders for saving
-#
-#y_model = np.array([])
-#w_Pt_vec = np.array([0.2, 0.1, 0.05, 0.025]) # Pt-loading values to loop through
-#for i, w in enumerate(w_Pt_vec):
-#    from pemfc_runner import *
-#    import pemfc_runner as user_inputs
-#    
-#    save = user_inputs.save = 0 # change to one if wanting to save outputs
-#    w_Pt = user_inputs.w_Pt = w
-#                
-#    folder_name = user_inputs.folder_name = folder_nm_generic + str(w_Pt) + 'mgcm-2'
-#    exec(open("Shared_Funcs/pemfc_pre.py").read())
-#    exec(open("Shared_Funcs/pemfc_dsvdt.py").read())
-#    exec(open("Shared_Funcs/pemfc_post.py").read())
-#    y_model = np.hstack([y_model, dphi_ss[1:]])
-#    if save == 1:
-#        ModuleWriter(cwd +'/Saved_Results/' +folder_name +'/user_inputs.csv', user_inputs)
+folder_nm_generic = 'cs_Pt_loop_' # pre-fix to all folders for saving
+
+y_model = np.array([])
+w_Pt_vec = np.array([0.2, 0.1, 0.05, 0.025]) # Pt-loading values to loop through
+for i, w in enumerate(w_Pt_vec):
+    from pemfc_runner import *
+    import pemfc_runner as user_inputs
+
+    save = user_inputs.save = 0 # change to one if wanting to save outputs
+    w_Pt = user_inputs.w_Pt = w
+
+    folder_name = user_inputs.folder_name = folder_nm_generic + str(w_Pt) + 'mgcm-2'
+    exec(open("Shared_Funcs/pemfc_pre.py").read())
+    exec(open("Shared_Funcs/pemfc_dsvdt.py").read())
+    exec(open("Shared_Funcs/pemfc_post.py").read())
+    y_model = np.hstack([y_model, dphi_ss[1:]])
+    if save == 1:
+        ModuleWriter(cwd +'/Saved_Results/' +folder_name +'/user_inputs.csv', user_inputs)
 
 " Plot sig_io and D_O2 as f(w_Pt) for 'mix' method "
 #clr = ['C0', 'C1', 'C2', 'C3']

Shared_Funcs/pemfc_post.py

@@ -38,7 +38,7 @@
 """ Post-processing for Plotting and Other Results """
 "-----------------------------------------------------------------------------"
 # Update plot settings:
-plt.close('all')
+#plt.close('all')
 font = plt.matplotlib.font_manager.FontProperties(family=font_nm, size=font_sz)
 plt.rcParams.update({'font.size': font_sz})
 
@@ -374,7 +374,7 @@
         yerr = np.array([0, 0.012, 0.007, 0.007, 0.012, 0.001, 0.008, 0.007,
                          0.007, 0.009, 0.009])
         color = 'C0'
-    elif all([data, w_Pt == 0.1}):
+    elif all([data, w_Pt == 0.1]):
         x = np.array([0.000, 0.050, 0.200, 0.400, 0.800, 1.000, 1.200, 1.500, 
                       1.650, 1.850, 2.000])
         y = np.array([0.930, 0.834, 0.785, 0.754, 0.711, 0.691, 0.673, 0.649, 

Shared_Funcs/pemfc_pre.py

@@ -16,7 +16,8 @@
 naf_b_ca = ct.Solution(ctifile, 'naf_bulk_ca')
 naf_b_ca.TP = T_ca, P_ca
 
-pt_s_ca = ct.Interface(ctifile, 'Pt_surf_ca_2s', [carb_ca, naf_b_ca, gas_ca])
+rm = str(rxn_mech)
+pt_s_ca = ct.Interface(ctifile, 'Pt_surf_ca_' + rm, [carb_ca, naf_b_ca, gas_ca])
 pt_s_ca.TP = T_ca, P_ca
 
 naf_s_ca = ct.Interface(ctifile, 'naf_surf_ca', [naf_b_ca, gas_ca])
@@ -52,8 +53,9 @@
 # Change parameters for optimization only:
 if 'optimize' in vars(): None
 else: optimize = 0
-
-if optimize == 1:
+
+# If running optimization w/ the 1-step reaction mechanism:
+if all([optimize == 1, rxn_mech == '1s']):
     if tog == 2:
         R_naf = (c*w_Pt**d + e) *1e-3
     else:
@@ -70,10 +72,31 @@
     pt_s_ca.modify_reaction(0, Rxn1)
 
     O2 = naf_b_ca.species(naf_b_ca.species_index('O2(Naf)'))
-    fuel_coeffs = O2.thermo.coeffs
-    fuel_coeffs[[1,8]] = np.array([3.28253784 +offset_opt, 3.782456360 +offset_opt])
+    O2_thermo = O2.thermo.coeffs
+    O2_thermo[[1,8]] = np.array([3.28253784 +offset_opt, 3.782456360 +offset_opt])
     O2.thermo = ct.NasaPoly2(200, 3500, ct.one_atm, fuel_coeffs)
     naf_b_ca.modify_species(naf_b_ca.species_index('O2(Naf)'), O2)
+
+# If running optimization w/ the 2-step reaction mechanism:
+elif all([optimize == 1, rxn_mech == '2s']):
+    R_naf = R_naf_opt
+    theta = theta_opt
+
+    A = A_opt
+    i_o = i_o_opt
+    Rxn1 = pt_s_ca.reaction(0)
+    Rxn2 = pt_s_ca.reaction(1)
+    Rxn1.rate = ct.Arrhenius(A, 0, 0)
+    Rxn2.rate = ct.Arrhenius(i_o, 0, 0)
+    pt_s_ca.modify_reaction(0, Rxn1)
+    pt_s_ca.modify_reaction(1, Rxn2)
+
+    O2 = naf_b_ca.species(naf_b_ca.species_index('O2(Naf)'))
+    O2_thermo = O2.thermo.coeffs
+    O2_thermo[[1,8]] = np.array([3.28253784 +offset_opt, 3.782456360 +offset_opt])
+    O2.thermo = ct.NasaPoly2(200, 3500, ct.one_atm, fuel_coeffs)
+    naf_b_ca.modify_species(naf_b_ca.species_index('O2(Naf)'), O2)
+
 ###############################################################################
 
 " Store phases in a common 'objs' dict: "

fitting_optimization.py

@@ -6,18 +6,19 @@
 Optimization routine:
     Minimize the error between the data from Owejan et. al. paper and the PEMFC
     models (core-shell and flooded-agglomerate) by varying the following:
-    theta, offset, R_naf, and i_o. 
+    theta, offset, R_naf, i_o, and A (if using a 2-step reaction mechanism). 
     
 Instructions:
     To use this script, open up the pemfc_runner.py file and comment out the 
     line of code that specifies w_Pt. Copy and paste the lines of code
-    below to input i_OCV, i_ext1, i_ext2, and i_ext3 into pemfc_runner.py. Once
-    these items have been completed, fill out the initial guess values for the 
-    optimization in 'p0', 'p1', or 'p2' in the order commented below. Ranges 
-    for each of these parameters should also be specified in the respective 
-    'bounds*' variable below. Simply run this script after these inputs have 
-    been filled out and allow for the error to be minimized for the data sets 
-    using the package scipy.optimize.minize.
+    below to input i_OCV, i_ext1, i_ext2, and i_ext3 into pemfc_runner.py. 
+    
+    Once these items have been completed, fill out the initial guess values for
+    the optimization in 'p*' in the order commented below. Ranges for each of 
+    these parameters should also be specified in the respective 'b*' variable 
+    below. Simply run this script after these inputs have been filled out and 
+    allow for the error to be minimized for the data sets using the package 
+    scipy.optimize.minize.
     
 Optimization options:
     The 'tog' variable can be set to 0, 1, or 2 in order to change which set of
@@ -32,9 +33,9 @@
     fitting parameters.
     
 Units for inputs:
-    R_naf [mOhm*cm^2], theta [degrees], i_o [A/cm^2 *(cm^3/kmol)^6.5], offset [-], 
-    a [A/cm^2 *(cm^3/kmol)^6.5 /(mg/cm^2)], b [A/cm^2 *(cm^3/kmol)^6.5], 
-    c [mOhm*cm^2 /(mg/cm^2)], and d [-], e [mOhm*cm^2].
+    R_naf [mOhm*cm^2], theta [degrees], i_o [A/cm^2 *(cm^3/kmol)^#], offset [-], 
+    a [A/cm^2 *(cm^3/kmol)^# /(mg/cm^2)], b [A/cm^2 *(cm^3/kmol)^#], 
+    c [mOhm*cm^2 /(mg/cm^2)], d [-], e [mOhm*cm^2], and A [(kmol/cm^3)^# / s].
 """
 
 import numpy as np
@@ -47,21 +48,31 @@
 i_ext2 = np.array([])
 i_ext3 = np.array([])"""
 
-# User Inputs for optimization of R_naf, theta, i_o, and offset:
-p0 = np.array([45e-3, 45, 4e-6, 0.55]) 
-bounds0 = np.array([[20e-3, 120e-3], [42.5, 80], [1e-6, 9e-6], [0.52, 1.0]]) # [min, max]
+###############################################################################
+""" Use this section if optimizing with a 1-step reaction mechanism. """
+
+# Inputs for optimization of R_naf, theta, i_o, and offset:
+p0_1s = np.array([45e-3, 45, 4e-6, 0.55]) 
+b0_1s = np.array([[20e-3, 120e-3], [42.5, 80], [1e-6, 9e-6], [0.52, 1.0]])
 
-# User Inputs for optimization of R_naf, theta, offset, a, and b from i_o = a*w_Pt + b:
-p1 = np.array([55e-3, 45, 0.55, 24.2e-6, -0.40e-6])
-bounds1 = np.array([[20e-3, 120e-3], [42.5, 80], [0.52, 1.0], [24.17e-6, 50e-6], [-0.5e-6, 1e-6]])
+# Inputs for optimization of R_naf, theta, offset, a, and b from i_o = a*w_Pt + b:
+p1_1s = np.array([55e-3, 45, 0.55, 24.2e-6, -0.40e-6])
+b1_1s = np.array([[20e-3, 120e-3], [42.5, 80], [0.52, 1.0], [24.17e-6, 50e-6], [-0.5e-6, 1e-6]])
 
-# User Inputs for optimization of theta, i_o, offset, c, d, e from R_naf = (c*w_Pt^d +e)*1e-3:
-p2 = np.array([50, 1.3e-6, 0.65, 0.929, -1.249, 35])
-bounds2 = np.array([[42.5, 80], [5e-7, 9e-6], [0.52, 1.0], [0., 1.5], [-1.5, 0.], [0., 120.]])
+# Inputs for optimization of theta, i_o, offset, c, d, and e from R_naf = (c*w_Pt^d +e)*1e-3:
+p2_1s = np.array([50, 1.3e-6, 0.65, 0.929, -1.249, 35])
+b2_1s = np.array([[42.5, 80], [5e-7, 9e-6], [0.52, 1.0], [0., 1.5], [-1.5, 0.], [0., 120.]])
 
 # Toggle for which optimization to run (0, 1, or 2):
 tog = 2
 
+###############################################################################
+""" Use this section if optimizing with a 2-step reaction mechanism. """
+
+# Inputs for optimization of R_naf, theta, i_o, offset, and A:
+p0_2s = np.array([])
+b0_2s = np.array([])
+
 """ Do not edit anything below this line """
 "-----------------------------------------------------------------------------"
 ###############################################################################
@@ -90,7 +101,7 @@
 y_data = np.hstack([y1[1:], y2[1:], y3[1:], y4[1:]])
 s_data = np.hstack([s1[1:], s2[1:], s3[1:], s4[1:]])
 
-def chi_sq_func0(p_opt): # p_opt[:] = R_naf, theta, i_o, offset
+def chi_sq_func0_1s(p_opt): # p_opt[:] = R_naf, theta, i_o, offset
 
     global w_Pt, R_naf_opt, theta_opt, i_o_opt, offset_opt
     R_naf_opt, theta_opt, i_o_opt, offset_opt = p_opt
@@ -126,7 +137,7 @@ def chi_sq_func0(p_opt): # p_opt[:] = R_naf, theta, i_o, offset
 i_o = 1.47659887e-6
 offset = 0.549173549"""
 
-def chi_sq_func1(p_opt): # p_opt[:] = R_naf, theta, offset, a, b
+def chi_sq_func1_1s(p_opt): # p_opt[:] = R_naf, theta, offset, a, b
     w_Pt_vec = np.array([0.2, 0.1, 0.05, 0.025])
 
     global w_Pt, R_naf_opt, theta_opt, offset_opt, a, b
@@ -165,7 +176,7 @@ def chi_sq_func1(p_opt): # p_opt[:] = R_naf, theta, offset, a, b
 a = 24.17e-6
 b = -0.434032963e-6"""
 
-def chi_sq_func2(p_opt): # p_opt[:] = theta, offset, i_o, c, d, e
+def chi_sq_func2_1s(p_opt): # p_opt[:] = theta, i_o, offset, c, d, e
     w_Pt_vec = np.array([0.2, 0.1, 0.05, 0.025])
 
     global w_Pt, theta_opt, i_o_opt, offset_opt, c, d, e
@@ -206,12 +217,34 @@ def chi_sq_func2(p_opt): # p_opt[:] = theta, offset, i_o, c, d, e
 d = -1.23061654
 e = 34.5429028"""
 
-if tog == 0:
-    res = minimize(chi_sq_func0, p0, method='L-BFGS-B', bounds=bounds0)
-elif tog == 1:
-    res = minimize(chi_sq_func1, p1, method='L-BFGS-B', bounds=bounds1)
-elif tog == 2:
-    res = minimize(chi_sq_func2, p2, method='L-BFGS-B', bounds=bounds2)
+def chi_sq_func0_2s(p_opt): # p_opt[:] = R_naf, theta, i_o, offset, A
+    w_Pt_vec = np.array([0.2, 0.1, 0.05, 0.025])
+
+    global w_Pt, R_naf_opt, theta_opt, i_o_opt, offset_opt, A_opt
+    theta_opt, i_o_opt, offset_opt, A_opt = p_opt
+    print('\n\nR_naf, theta, i_o, offset, A:', p_opt)
+
+    y_model = np.array([])
+    for i, w in enumerate(w_Pt_vec):
+        w_Pt = w
+        exec(open("pemfc_runner.py").read(), globals(), globals())
+        y_model = np.hstack([y_model, dphi_ss[1:]])
+
+    chi_sq = sum((y_data - y_model)**2 / s_data**2)
+    print('\nchi_sq:', chi_sq)
+
+    return chi_sq
+
+import pemfc_runner as run
+
+if all([rxn_mech == '1s', tog == 0]):
+    res = minimize(chi_sq_func0_1s, p0_1s, method='L-BFGS-B', bounds=b0_1s)
+elif all([rxn_mech == '1s', tog == 1]):
+    res = minimize(chi_sq_func1_1s, p1_1s, method='L-BFGS-B', bounds=b1_1s)
+elif all([rxn_mech == '1s', tog == 2]):
+    res = minimize(chi_sq_func2_1s, p2_1s, method='L-BFGS-B', bounds=b2_1s)
+elif all([rxn_mech == '2s']):
+    res = minimize(chi_sq_func0_2s, p0_2s, method='L-BFGS-B', bounds=b0_2s)
 
 # Use this in w_Pt_runner to see what plots look like:
 """ Post this line of code in w_Pt_runner.py between line 38 and line 40. Set

pemfc_runner.py

@@ -71,9 +71,8 @@
 
 """ User Input Parameters """
 "-----------------------------------------------------------------------------"
-model = 'flooded_agg'               # CL geom: 'core_shell' or 'flooded_agg'
-ctifile = 'pemfc_cs.cti'            # cantera input file to match chosen model
-ver = 1                             # debugging radial diffusion (1:cs, 2:fa)
+model = 'core_shell'                # CL geom: 'core_shell' or 'flooded_agg'
+rxn_mech = '2s'                     # '1s' or '2s' for 1- or 2-step reaction
 
 " Initial electrochemical values "
 i_OCV = 0                           # 0 [A/cm^2] -> or single run if != 0 
@@ -94,12 +93,12 @@
 C_dl = 1.5e-9        # capacitance of double layer [F/m^2]
 sig_method = 'mix'   # 'lam', 'bulk', 'mix', or 'sun' for conductivity method
 D_O2_method = 'mix'  # 'lam', 'bulk', 'mix', or 'sun' for D_O2 method
-R_naf = 45e-3        # resistance of Nafion membrane [Ohm*cm^2] (45:cs, 60:fa)
+R_naf = 45e-3        # resistance of Nafion membrane [Ohm*cm^2]
 
 " Pt loading and geometric values "
 area_calcs = 1      # control for area calculations [0:p_Pt, 1:Pt_loading]
 p_Pt = 10           # percentage of Pt covering the carbon particle surface [%]
-#w_Pt = 0.2          # loading of Pt on carbon [mg/cm^2]
+w_Pt = 0.2          # loading of Pt on carbon [mg/cm^2]
 rho_Pt = 21.45e3    # density of Pt for use in area property calcs [kg/m^3]
 r_c = 50e-9         # radius of single carbon particle [m] (50:cs, 25:fa)
 r_Pt = 1e-9         # radius of Pt 1/2 sphere sitting on C surface [m]
@@ -144,11 +143,10 @@
 polar = 1           # turn on to generate full cell polarization curves
 over_p = 0          # turn on to plot overpotential curve for cathode
 
-<<<<<<< HEAD
-" Verification settings - (0: off and 1: on) unless otherwise stated "
-data = 1            # include data from Owejan et. al. on polarization if available
-i_ver = 1           # verify current between GDL and CL with O2 flux calcs
-i_find = 0.5        # current from polarization curve to use in i_ver processing
+" Verification settings "
+data = 1            # include data from Owejan et. al. if available (0: off, 1: on)
+i_ver = 1           # verify current between GDL and CL (0: off, 1: on)
+i_find = 0.5        # current from polarization to use in i_ver calculation [A/cm^2]
 
 " Plotting options "
 font_nm = 'Arial'   # name of font for plots