Progress toward heating bug fix

ecm/__funcs__.py

@@ -50,6 +50,7 @@ def adjust_parameters(parameter_values, I_typical):
             # "Typical current [A]": I_typical,
             "Current function [A]": current_function,
             "Electrode height [m]": "[input]",
+            # "Cell volume [m3]": "[input]",
         }
     )
     parameter_values.update({"Current": "[input]"}, check_already_exists=False)
@@ -76,6 +77,7 @@ def output_variables():
         "Volume-averaged irreversible electrochemical heating [W.m-3]",
         "X-averaged irreversible electrochemical heating [W.m-3]",
         "X-averaged reversible heating [W.m-3]",
+        "X-averaged total heating [W.m-3]",
         "Volume-averaged reversible heating [W.m-3]",
         'X-averaged Ohmic heating [W.m-3]',
         'Volume-averaged Ohmic heating [W.m-3]',
@@ -85,7 +87,11 @@ def output_variables():
         'Volume-averaged cell temperature [C]',
         'X-averaged negative electrode temperature [C]',
         'X-averaged separator temperature [C]',
-        'X-averaged positive electrode temperature [C]'
+        'X-averaged positive electrode temperature [C]',
+        "Total heating [W]",
+        "Total heating per unit electrode-pair area [W.m-2]",
+        "Volume-averaged total heating [W.m-3]",
+        "X-averaged total heating [W.m-3]",
     ]
 
 
@@ -242,22 +248,30 @@ def apply_heat_source_lp(project, Q):
     # The SPMs are defined at the throat but the pores represent the
     # Actual electrode volume so need to interpolate for heat sources
     phys = project.physics()["phys_01"]
-    phys["throat.heat_source"] = Q
+    phys["throat.heat_source"] = Q * 2
     phys.add_model(
         propname="pore.heat_source",
         model=op.models.misc.from_neighbor_throats,
         prop="throat.heat_source",
         mode="mean",
     )
 
+    # phys['pore.heat_source'] = phys.interpolate_data(propname='throat.heat_source')
+
+    # print('pore heat', phys['pore.heat_source'].sum())
+    # print('throat heat', phys['throat.heat_source'].sum())
+
 
 def run_step_transient(project, time_step, BC_value, cp, rho, third=False, **kwargs):
     # To Do - test whether this needs to be transient
     net = project.network
     phase = project.phases()["phase_01"]
     phys = project.physics()["phys_01"]
     phys["pore.A1"] = 0.0
-    Q_spm = phys["pore.heat_source"] * net["pore.volume"]
+    Q_spm = phys["pore.heat_source"] * net["pore.volume"] * len(net['pore.volume'])
+    # print(Q_spm.max(), Q_spm.min())
+    # print('total W heat produced',Q_spm.sum())
+    # print('total pore volume', net["pore.volume"].sum())
     # Q_cc = net["pore.cc_power_loss"]
     # print(
     #     "Q_spm",

ecm/__liionsolve__.py

@@ -38,8 +38,8 @@ def do_heating():
 
 def run_simulation_lp(parameter_values, experiment, initial_soc, project, **kwargs):
     ###########################################################################
-    kwargs.setdefault('t_slice', 10)
-    kwargs.setdefault('t_precision', 12)
+    kwargs.setdefault('t_slice', 4)
+    kwargs.setdefault('t_precision', 1)
     kwargs.setdefault('disable', True)
     kwargs.setdefault('max_workers', 1)
     ###########################################################################
@@ -131,6 +131,7 @@ def run_simulation_lp(parameter_values, experiment, initial_soc, project, **kwar
         "Electrode height [m]": e_heights,
         "Input temperature [K]": spm_temperature
     }
+    # print('min e height', e_heights.min(), 'max e height', e_heights.max())
     ###########################################################################
     # Initialisation
     experiment_init = pybamm.Experiment(
@@ -152,7 +153,8 @@ def run_simulation_lp(parameter_values, experiment, initial_soc, project, **kwar
         initial_soc=initial_soc,
         setup_only=True,
     )
-    Qvar = "Volume-averaged total heating [W.m-3]"
+    # Qvar = "Volume-averaged total heating [W.m-3]"
+    Qvar = "X-averaged total heating [W.m-3]"
     Qid = np.argwhere(np.asarray(manager.variable_names) == Qvar).flatten()[0]
     lp.logger.notice("Starting initial step solve")
     vlims_ok = True
@@ -169,7 +171,9 @@ def run_simulation_lp(parameter_values, experiment, initial_soc, project, **kwar
             ###################################################################
             # Apply Heat Sources
             Q_tot = manager.output[Qid, step, :]
+            # print('Qtot sum', Q_tot.sum())
             Q = get_cc_power_loss(net, netlist)
+            # print('cc power loss sum', Q.sum())
             # print(Q_tot)
             # print(Q)
             # To do - Get cc heat from netlist
@@ -185,6 +189,7 @@ def run_simulation_lp(parameter_values, experiment, initial_soc, project, **kwar
                 project, dim_time_step, T0, cp, rho, thermal_third, **kwargs
             )
             # Interpolate the node temperatures for the SPMs
+            print(f'before: (Max)[{np.round(spm_temperature.max(),2)}, (Min)[{np.round(spm_temperature.min(),2)}]]')
             spm_temperature = phase.interpolate_data("pore.temperature")[res_Ts]
             # T_non_dim_spm = fT_non_dim(parameter_values, spm_temperature)
             ###################################################################
@@ -218,7 +223,8 @@ def run_simulation_lp(parameter_values, experiment, initial_soc, project, **kwar
         initial_soc=initial_soc,
         setup_only=True,
     )
-    Qvar = "Volume-averaged total heating [W.m-3]"
+    # Qvar = "Volume-averaged total heating [W.m-3]"
+    Qvar = "X-averaged total heating [W.m-3]"
     Qid = np.argwhere(np.asarray(manager.variable_names) == Qvar).flatten()[0]
     lp.logger.notice("Starting step solve")
     vlims_ok = True
@@ -234,19 +240,25 @@ def run_simulation_lp(parameter_values, experiment, initial_soc, project, **kwar
             ###################################################################
             # Apply Heat Sources
             Q_tot = manager.output[Qid, step, :]
+            # print('Qtot sum', Q_tot.sum())
             Q = get_cc_power_loss(net, netlist)
+            # print('cc power loss sum', Q.sum())
             # To do - Get cc heat from netlist
             # Q_ohm_cc = net.interpolate_data("pore.cc_power_loss")[res_Ts]
             # Q_ohm_cc /= net["throat.volume"][res_Ts]
             # key = "Volume-averaged Ohmic heating CC [W.m-3]"
             # vh[key][outer_step, :] = Q_ohm_cc[sorted_res_Ts]
             Q[res_Ts] += Q_tot
+            # print('Q mean', Q.mean())
+            # print('volume sum', net["pore.volume"].sum())
+            print('total W heat produced pre-solve', Q.sum() * net["pore.volume"].sum())
             ecm.apply_heat_source_lp(project, Q)
             # Calculate Global Temperature
             ecm.run_step_transient(
                 project, dim_time_step, T0, cp, rho, thermal_third, **kwargs
             )
             # Interpolate the node temperatures for the SPMs
+            print(f'before: (Max)[{np.round(spm_temperature.max(),5)}, (Min)[{np.round(spm_temperature.min(),5)}]')
             spm_temperature = phase.interpolate_data("pore.temperature")[res_Ts]
             ###################################################################
             if vlims_ok:

ecm/autofit_jellyroll.py

@@ -135,6 +135,8 @@ def get_spiral_params(parameter_values, form_factor='18650', positive_tab=None,
     else:
         raise ValueError("form_factor must be '18650', '21700' or 'pouch'")
 
+    # parameter_values.update({'Electrode height [m]': length3d})
+
     Nlayers, L = calculate_spiral(inner_diameter, outer_diameter, spacing)
 
     # nominal_area = parameter_values['Electrode height [m]'] * parameter_values['Electrode width [m]']

ecm/runfuns.py

@@ -73,6 +73,8 @@ def run_simulation_cylindrical(i, row):
     else:
         raise ValueError("form_factor must be '18650' or 'pouch'")
 
+    # parameter_values.update({'Electrode height [m]': 0.00065})
+    # parameter_values.update({'Electrode width [m]': 0.065})
     project, arc_edges = get_spiral_params(parameter_values)
 
     # ecm.plot_topology(project.network)
@@ -86,16 +88,15 @@ def run_simulation_cylindrical(i, row):
     # get the experiment
     experiment = get_experiment(row,estimated_capacity)
 
-    parameter_values.update({'Electrode width [m]': 0.065})
 
-    trans_kwargs = {'t_slice':10,'t_precision':1}
+    # trans_kwargs = {'t_slice':10,'t_precision':1}
 
     # Run simulation
     project, output = ecm.run_simulation_lp(parameter_values=parameter_values,
                                             experiment=experiment,
                                             initial_soc=None,
                                             project=project,
-                                            **trans_kwargs)
+                                            )
 
 
 
@@ -166,7 +167,7 @@ def get_experiment(row,capacity):
             [
                 ('Rest for 1 second'),
                 (f"Discharge at {row['rate']*capacity} A for 10 seconds")
-            ], period='0.5 seconds'
+            ], period='0.05 seconds'
         )
     elif row['experiment_type'] == 'SOC DCR':
         experiment = pybamm.Experiment(
@@ -176,10 +177,11 @@ def get_experiment(row,capacity):
             ], period='0.5 seconds'
         )
     elif row['experiment_type'] == 'Full Discharge':
-        steps = 1/row['rate']*3600/360
+        seconds_Iapp = np.round(1/row['rate']*3600,5)
+        steps = seconds_Iapp/360
         experiment = pybamm.Experiment(
             [
-                (f"Discharge at {row['rate']*capacity} A for {1/row['rate']*3600} seconds")
+                (f"Discharge at {row['rate']*np.round(capacity,5)} A for {seconds_Iapp} seconds")
             ], period=f'{steps} seconds'
         )
     elif row['experiment_type'] == 'Discharge DC':
@@ -226,6 +228,12 @@ def get_experiment(row,capacity):
                 (f"Discharge at {0.2*capacity} A for 30 minutes (10 s period)", f"Discharge at {2*capacity} A for 30 seconds or until 2.5V (1 s period)", f"Charge at {0.2*capacity} A for 5 minutes (1 s period)"),
             ]
         )
+    elif row['experiment_type'] == 'Rest 1 hour':
+        experiment = pybamm.Experiment(
+            [
+                (f"Rest for 1 hour (30 second period)"),
+            ]
+        )
     else:
         raise ValueError(f"Unknown experiment type: {row['experiment_type']}")
     return experiment