From 5143e07f9c82e9e40fb680bd1e6b807327c5a5e2 Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Fri, 25 Oct 2019 09:04:16 +0100 Subject: [PATCH 01/11] #678 fix dimension typo --- .../submodels/electrolyte/base_electrolyte_conductivity.py | 2 +- pybamm/parameters/standard_parameters_lead_acid.py | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py b/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py index 7d685dc7fe..cb3c49e64f 100644 --- a/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py +++ b/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py @@ -250,7 +250,7 @@ def _get_domain_current_variables(self, i_e, domain=None): variables = { domain + " electrolyte current density": i_e, - domain + " electrolyte current density [V]": i_e * i_typ, + domain + " electrolyte current density [A.m-2]": i_e * i_typ, } return variables diff --git a/pybamm/parameters/standard_parameters_lead_acid.py b/pybamm/parameters/standard_parameters_lead_acid.py index 90840bd316..496d5fb85e 100644 --- a/pybamm/parameters/standard_parameters_lead_acid.py +++ b/pybamm/parameters/standard_parameters_lead_acid.py @@ -416,6 +416,9 @@ def U_p_dimensional(c_e, T): c_n_init = c_e_init c_p_init = c_e_init +# Thermal effects not implemented for lead-acid, but parameter needed for consistency +Theta = pybamm.Scalar(0) # ratio of typical temperature change to ambient temperature + # -------------------------------------------------------------------------------------- "5. Dimensionless Functions" From 84206d59a13fc0b56ef051aa3f7e7e7e029a11e6 Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Mon, 28 Oct 2019 14:06:28 +0000 Subject: [PATCH 02/11] #678 thermal comparison --- .../compare_comsol_thermal.py | 126 ++++++++++++++++ results/comsol_comparison/load_comsol_data.py | 138 ++++++++++++++++++ 2 files changed, 264 insertions(+) create mode 100644 results/comsol_comparison/compare_comsol_thermal.py create mode 100644 results/comsol_comparison/load_comsol_data.py diff --git a/results/comsol_comparison/compare_comsol_thermal.py b/results/comsol_comparison/compare_comsol_thermal.py new file mode 100644 index 0000000000..7c8d5d49a2 --- /dev/null +++ b/results/comsol_comparison/compare_comsol_thermal.py @@ -0,0 +1,126 @@ +import pybamm +import numpy as np +import os +import pickle +import scipy.interpolate as interp +import matplotlib.pyplot as plt + +# change working directory to the root of pybamm +os.chdir(pybamm.root_dir()) + +"-----------------------------------------------------------------------------" +"Load comsol data" + +comsol_variables = pickle.load( + open("input/comsol_results/comsol_thermal_1C.pickle", "rb") +) + +"-----------------------------------------------------------------------------" +"Create and solve pybamm model" + +# load model and geometry +pybamm.set_logging_level("INFO") +options = {"thermal": "x-full"} +pybamm_model = pybamm.lithium_ion.DFN(options) +geometry = pybamm_model.default_geometry + +# load parameters and process model and geometry +param = pybamm_model.default_parameter_values +param.process_model(pybamm_model) +param.process_geometry(geometry) + +# create mesh +var = pybamm.standard_spatial_vars +var_pts = {var.x_n: 31, var.x_s: 11, var.x_p: 31, var.r_n: 11, var.r_p: 11} +mesh = pybamm.Mesh(geometry, pybamm_model.default_submesh_types, var_pts) + +# discretise model +disc = pybamm.Discretisation(mesh, pybamm_model.default_spatial_methods) +disc.process_model(pybamm_model) + +# discharge timescale +tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge) + +# solve model at comsol times +time = comsol_variables["time"] / tau.evaluate(0) +solution = pybamm_model.default_solver.solve(pybamm_model, time) + +"-----------------------------------------------------------------------------" +"Make Comsol 'model' for comparison" + +whole_cell = ["negative electrode", "separator", "positive electrode"] +comsol_t = comsol_variables["time"] +L_x = param.evaluate(pybamm.standard_parameters_lithium_ion.L_x) + + +def get_interp_fun(variable, domain): + """ + Create a :class:`pybamm.Function` object using the variable, to allow plotting with + :class:`'pybamm.QuickPlot'` (interpolate in space to match edges, and then create + function to interpolate in time) + """ + if domain == ["negative electrode"]: + comsol_x = comsol_variables["x_n"] + elif domain == ["separator"]: + comsol_x = comsol_variables["x_s"] + elif domain == ["positive electrode"]: + comsol_x = comsol_variables["x_p"] + elif domain == whole_cell: + comsol_x = comsol_variables["x"] + # Make sure to use dimensional space + pybamm_x = mesh.combine_submeshes(*domain)[0].nodes * L_x + variable = interp.interp1d(comsol_x, variable, axis=0)(pybamm_x) + + def myinterp(t): + return interp.interp1d(comsol_t, variable)(t)[:, np.newaxis] + + # Make sure to use dimensional time + fun = pybamm.Function(myinterp, pybamm.t * tau) + fun.domain = domain + return fun + + +comsol_c_n_surf = get_interp_fun(comsol_variables["c_n_surf"], ["negative electrode"]) +comsol_c_e = get_interp_fun(comsol_variables["c_e"], whole_cell) +comsol_c_p_surf = get_interp_fun(comsol_variables["c_p_surf"], ["positive electrode"]) +comsol_phi_n = get_interp_fun(comsol_variables["phi_n"], ["negative electrode"]) +comsol_phi_e = get_interp_fun(comsol_variables["phi_e"], whole_cell) +comsol_phi_p = get_interp_fun(comsol_variables["phi_p"], ["positive electrode"]) +comsol_voltage = interp.interp1d(comsol_t, comsol_variables["voltage"]) +comsol_temperature = get_interp_fun(comsol_variables["temperature"], whole_cell) +comsol_q_irrev_n = get_interp_fun(comsol_variables["Q_irrev_n"], ["negative electrode"]) +comsol_q_irrev_p = get_interp_fun(comsol_variables["Q_irrev_p"], ["positive electrode"]) +comsol_q_rev_n = get_interp_fun(comsol_variables["Q_rev_n"], ["negative electrode"]) +comsol_q_rev_p = get_interp_fun(comsol_variables["Q_rev_p"], ["positive electrode"]) +comsol_q_total_n = get_interp_fun(comsol_variables["Q_total_n"], ["negative electrode"]) +comsol_q_total_s = get_interp_fun(comsol_variables["Q_total_s"], ["separator"]) +comsol_q_total_p = get_interp_fun(comsol_variables["Q_total_p"], ["positive electrode"]) + +# Create comsol model with dictionary of Matrix variables +comsol_model = pybamm.BaseModel() +comsol_model.variables = { + "Negative particle surface concentration [mol.m-3]": comsol_c_n_surf, + "Electrolyte concentration [mol.m-3]": comsol_c_e, + "Positive particle surface concentration [mol.m-3]": comsol_c_p_surf, + "Current [A]": pybamm_model.variables["Current [A]"], + "Negative electrode potential [V]": comsol_phi_n, + "Electrolyte potential [V]": comsol_phi_e, + "Positive electrode potential [V]": comsol_phi_p, + "Terminal voltage [V]": pybamm.Function(comsol_voltage, pybamm.t * tau), + "Cell temperature [K]": comsol_temperature, +} + +"-----------------------------------------------------------------------------" +"Plot comparison" + +#TODO: fix QuickPlot +for var in comsol_model.variables.keys(): + plot = pybamm.QuickPlot( + [pybamm_model, comsol_model], + mesh, + [solution, solution], + output_variables=[var], + labels=["PyBaMM", "Comsol"], + ) + plot.dynamic_plot(testing=True) +plt.show() diff --git a/results/comsol_comparison/load_comsol_data.py b/results/comsol_comparison/load_comsol_data.py new file mode 100644 index 0000000000..ea6f8366e0 --- /dev/null +++ b/results/comsol_comparison/load_comsol_data.py @@ -0,0 +1,138 @@ +import pybamm +import os +import pandas as pd +import pickle +import numpy as np + +# change working directory the root of pybamm +os.chdir(pybamm.root_dir()) + +# set filepath for data and name of file to pickle to +path = "input/comsol_results_csv/thermal/1C/" +savefile = "input/comsol_results/comsol_thermal_1C.pickle" + +# time-voltage +comsol = pd.read_csv(path + "voltage.csv", sep=",", header=None) +comsol_time = comsol[0].values +comsol_time_npts = len(comsol_time) +comsol_voltage = comsol[1].values + +# negative electrode potential +comsol = pd.read_csv(path + "phi_s_n.csv", sep=",", header=None) +comsol_x_n_npts = int(len(comsol[0].values) / comsol_time_npts) +comsol_x_n = comsol[0].values[0:comsol_x_n_npts] +comsol_phi_n_vals = np.reshape( + comsol[1].values, (comsol_x_n_npts, comsol_time_npts), order="F" +) + +# negative particle surface concentration +comsol = pd.read_csv(path + "c_n_surf.csv", sep=",", header=None) +comsol_c_n_surf_vals = np.reshape( + comsol[1].values, (comsol_x_n_npts, comsol_time_npts), order="F" +) + +# positive electrode potential +comsol = pd.read_csv(path + "phi_s_p.csv", sep=",", header=None) +comsol_x_p_npts = int(len(comsol[0].values) / comsol_time_npts) +comsol_x_p = comsol[0].values[0:comsol_x_p_npts] +comsol_phi_p_vals = np.reshape( + comsol[1].values, (comsol_x_p_npts, comsol_time_npts), order="F" +) + +# positive particle surface concentration +comsol = pd.read_csv(path + "c_p_surf.csv", sep=",", header=None) +comsol_c_p_surf_vals = np.reshape( + comsol[1].values, (comsol_x_p_npts, comsol_time_npts), order="F" +) + +# electrolyte concentration +comsol = pd.read_csv(path + "c_e.csv", sep=",", header=None) +comsol_x_npts = int(len(comsol[0].values) / comsol_time_npts) +comsol_x = comsol[0].values[0:comsol_x_npts] +comsol_c_e_vals = np.reshape( + comsol[1].values, (comsol_x_npts, comsol_time_npts), order="F" +) + +# electrolyte potential +comsol = pd.read_csv(path + "phi_e.csv", sep=",", header=None) +comsol_phi_e_vals = np.reshape( + comsol[1].values, (comsol_x_npts, comsol_time_npts), order="F" +) + +# temperature +comsol = pd.read_csv(path + "temperature.csv", sep=",", header=None) +comsol_temp_vals = np.reshape( + comsol[1].values, (comsol_x_npts, comsol_time_npts), order="F" +) + +# irreversible heat source in negative electrode +comsol = pd.read_csv(path + "q_irrev_n.csv", sep=",", header=None) +comsol_q_irrev_n_vals = np.reshape( + comsol[1].values, (comsol_x_n_npts, comsol_time_npts), order="F" +) + +# irreversible heat source in positive electrode +comsol = pd.read_csv(path + "q_irrev_p.csv", sep=",", header=None) +comsol_q_irrev_p_vals = np.reshape( + comsol[1].values, (comsol_x_p_npts, comsol_time_npts), order="F" +) + +# reversible heat source in negative electrode +comsol = pd.read_csv(path + "q_rev_n.csv", sep=",", header=None) +comsol_q_rev_n_vals = np.reshape( + comsol[1].values, (comsol_x_n_npts, comsol_time_npts), order="F" +) + +# reversible heat source in positive electrode +comsol = pd.read_csv(path + "q_rev_p.csv", sep=",", header=None) +comsol_q_rev_p_vals = np.reshape( + comsol[1].values, (comsol_x_p_npts, comsol_time_npts), order="F" +) + +# total heat source in negative electrode +comsol = pd.read_csv(path + "q_total_n.csv", sep=",", header=None) +comsol_q_total_n_vals = np.reshape( + comsol[1].values, (comsol_x_n_npts, comsol_time_npts), order="F" +) + +# total heat source in separator +comsol = pd.read_csv(path + "q_total_s.csv", sep=",", header=None) +comsol_x_s_npts = int(len(comsol[0].values) / comsol_time_npts) +comsol_x_s = comsol[0].values[0:comsol_x_s_npts] +comsol_q_total_s_vals = np.reshape( + comsol[1].values, (comsol_x_s_npts, comsol_time_npts), order="F" +) + +# total heat source in positive electrode +comsol = pd.read_csv(path + "q_total_p.csv", sep=",", header=None) +comsol_q_total_p_vals = np.reshape( + comsol[1].values, (comsol_x_p_npts, comsol_time_npts), order="F" +) + + +# add comsol variables to dict and pickle +comsol_variables = { + "time": comsol_time, + "x_n": comsol_x_n, + "x_s": comsol_x_s, + "x_p": comsol_x_p, + "x": comsol_x, + "voltage": comsol_voltage, + "phi_n": comsol_phi_n_vals, + "phi_p": comsol_phi_p_vals, + "phi_e": comsol_phi_e_vals, + "c_n_surf": comsol_c_n_surf_vals, + "c_p_surf": comsol_c_p_surf_vals, + "c_e": comsol_c_e_vals, + "temperature": comsol_temp_vals, + "Q_irrev_n": comsol_q_irrev_n_vals, + "Q_irrev_p": comsol_q_irrev_p_vals, + "Q_rev_n": comsol_q_rev_n_vals, + "Q_rev_p": comsol_q_rev_p_vals, + "Q_total_n": comsol_q_total_n_vals, + "Q_total_s": comsol_q_total_s_vals, + "Q_total_p": comsol_q_total_p_vals, +} + +with open(savefile, "wb") as f: + pickle.dump(comsol_variables, f, pickle.HIGHEST_PROTOCOL) From ca01dd7dabb9d717e4ab90d3394b0dc97007fff6 Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Tue, 29 Oct 2019 17:23:49 +0000 Subject: [PATCH 03/11] #678 fix cooling --- examples/notebooks/models/SPM.ipynb | 12 +- examples/scripts/thermal_lithium_ion.py | 5 +- .../models/submodels/thermal/base_thermal.py | 249 ++++++++++++++---- .../thermal/isothermal/isothermal.py | 12 +- .../x_full/x_full_no_current_collector.py | 6 +- .../thermal/x_lumped/base_x_lumped.py | 8 + .../x_lumped_0D_current_collectors.py | 8 +- .../x_lumped_1D_current_collectors.py | 9 +- .../x_lumped_2D_current_collectors.py | 5 +- .../x_lumped_no_current_collectors.py | 46 +++- .../xyz_lumped_1D_current_collector.py | 1 + .../xyz_lumped_2D_current_collector.py | 1 + pybamm/parameters/geometric_parameters.py | 4 +- .../standard_parameters_lithium_ion.py | 3 +- pybamm/parameters/thermal_parameters.py | 16 +- .../compare_comsol_thermal.py | 215 ++++++++++++++- results/comsol_comparison/load_comsol_data.py | 4 + 17 files changed, 499 insertions(+), 105 deletions(-) diff --git a/examples/notebooks/models/SPM.ipynb b/examples/notebooks/models/SPM.ipynb index 54bf98ee01..47a3dd06c0 100644 --- a/examples/notebooks/models/SPM.ipynb +++ b/examples/notebooks/models/SPM.ipynb @@ -638,17 +638,17 @@ "\t- Local current collector potential difference\n", "\t- Local current collector potential difference [V]\n", "\t- Ohmic heating\n", - "\t- Ohmic heating [A.V.m-3]\n", + "\t- Ohmic heating [W.m-3]\n", "\t- Irreversible electrochemical heating\n", - "\t- Irreversible electrochemical heating [A.V.m-3]\n", + "\t- Irreversible electrochemical heating [W.m-3]\n", "\t- Reversible heating\n", - "\t- Reversible heating [A.V.m-3]\n", + "\t- Reversible heating [W.m-3]\n", "\t- Total heating\n", - "\t- Total heating [A.V.m-3]\n", + "\t- Total heating [W.m-3]\n", "\t- X-averaged total heating\n", - "\t- X-averaged total heating [A.V.m-3]\n", + "\t- X-averaged total heating [W.m-3]\n", "\t- Volume-averaged total heating\n", - "\t- Volume-averaged total heating [A.V.m-3]\n", + "\t- Volume-averaged total heating [W.m-3]\n", "\t- Positive current collector potential [V]\n", "\t- Local potential difference\n", "\t- Local potential difference [V]\n", diff --git a/examples/scripts/thermal_lithium_ion.py b/examples/scripts/thermal_lithium_ion.py index 13f3826eb2..79747040d8 100644 --- a/examples/scripts/thermal_lithium_ion.py +++ b/examples/scripts/thermal_lithium_ion.py @@ -18,7 +18,7 @@ # load parameter values and process models and geometry param = models[0].default_parameter_values -param.update({"Heat transfer coefficient [W.m-2.K-1]": 0.1}) +param.update({"Heat transfer coefficient [W.m-2.K-1]": 1}) for model in models: param.process_model(model) @@ -40,7 +40,8 @@ solutions = [None] * len(models) t_eval = np.linspace(0, 0.17, 100) for i, model in enumerate(models): - solution = model.default_solver.solve(model, t_eval) + solver = pybamm.ScipySolver(atol=1e-8, rtol=1e-8) + solution = solver.solve(model, t_eval) solutions[i] = solution # plot diff --git a/pybamm/models/submodels/thermal/base_thermal.py b/pybamm/models/submodels/thermal/base_thermal.py index e2b801db14..4495694c82 100644 --- a/pybamm/models/submodels/thermal/base_thermal.py +++ b/pybamm/models/submodels/thermal/base_thermal.py @@ -92,16 +92,31 @@ def _get_standard_coupled_variables(self, variables): phi_s_n = variables["Negative electrode potential"] phi_s_p = variables["Positive electrode potential"] + # Ohmic heating in solid Q_ohm_s_cn, Q_ohm_s_cp = self._current_collector_heating(variables) Q_ohm_s_n = -pybamm.inner(i_s_n, pybamm.grad(phi_s_n)) Q_ohm_s_s = pybamm.FullBroadcast(0, ["separator"], "current collector") Q_ohm_s_p = -pybamm.inner(i_s_p, pybamm.grad(phi_s_p)) Q_ohm_s = pybamm.Concatenation(Q_ohm_s_n, Q_ohm_s_s, Q_ohm_s_p) - Q_ohm_e = -pybamm.inner(i_e, pybamm.grad(phi_e)) + # Ohmic heating in electrolyte + # TODO: change full stefan-maxwell conductivity so that i_e is always + # a Concatenation + if isinstance(i_e, pybamm.Concatenation): + # compute by domain if possible + i_e_n, i_e_s, i_e_p = i_e.orphans + phi_e_n, phi_e_s, phi_e_p = phi_e.orphans + Q_ohm_e_n = -pybamm.inner(i_e_n, pybamm.grad(phi_e_n)) + Q_ohm_e_s = -pybamm.inner(i_e_s, pybamm.grad(phi_e_s)) + Q_ohm_e_p = -pybamm.inner(i_e_p, pybamm.grad(phi_e_p)) + Q_ohm_e = pybamm.Concatenation(Q_ohm_e_n, Q_ohm_e_s, Q_ohm_e_p) + else: + Q_ohm_e = -pybamm.inner(i_e, pybamm.grad(phi_e)) + # Total Ohmic heating Q_ohm = Q_ohm_s + Q_ohm_e + # Irreversible electrochemical heating Q_rxn_n = j_n * eta_r_n Q_rxn_p = j_p * eta_r_p Q_rxn = pybamm.Concatenation( @@ -112,6 +127,7 @@ def _get_standard_coupled_variables(self, variables): ] ) + # Reversible electrochemical heating Q_rev_n = j_n * (param.Theta ** (-1) + T_n) * dUdT_n Q_rev_p = j_p * (param.Theta ** (-1) + T_p) * dUdT_p Q_rev = pybamm.Concatenation( @@ -122,6 +138,7 @@ def _get_standard_coupled_variables(self, variables): ] ) + # Total heating Q = Q_ohm + Q_rxn + Q_rev # Compute the X-average over the current collectors by default. @@ -133,77 +150,194 @@ def _get_standard_coupled_variables(self, variables): Q_rev_av = self._x_average(Q_rev, 0, 0) Q_av = self._x_average(Q, Q_ohm_s_cn, Q_ohm_s_cp) + # Compute volume-averaged heat source terms Q_ohm_vol_av = self._yz_average(Q_ohm_av) Q_rxn_vol_av = self._yz_average(Q_rxn_av) Q_rev_vol_av = self._yz_average(Q_rev_av) Q_vol_av = self._yz_average(Q_av) + # Dimensional scaling for heat source terms + Q_scale = param.i_typ * param.potential_scale / param.L_x + + # Geat heat source by domain + if isinstance(Q_ohm_e, pybamm.Concatenation): + variables.update( + self._get_domain_heat_source_terms(Q_ohm_s, Q_ohm_e, Q_rxn, Q_rev) + ) + variables.update( { "Ohmic heating": Q_ohm, - "Ohmic heating [A.V.m-3]": param.i_typ - * param.potential_scale - * Q_ohm - / param.L_x, + "Ohmic heating [W.m-3]": Q_ohm * Q_scale, "X-averaged Ohmic heating": Q_ohm_av, - "X-averaged Ohmic heating [A.V.m-3]": param.i_typ - * param.potential_scale - * Q_ohm_av - / param.L_x, + "X-averaged Ohmic heating [W.m-3]": Q_ohm_av * Q_scale, "Volume-averaged Ohmic heating": Q_ohm_vol_av, - "Volume-averaged Ohmic heating [A.V.m-3]": param.i_typ - * param.potential_scale - * Q_ohm_vol_av - / param.L_x, + "Volume-averaged Ohmic heating [W.m-3]": Q_ohm_vol_av * Q_scale, "Irreversible electrochemical heating": Q_rxn, - "Irreversible electrochemical heating [A.V.m-3]": param.i_typ - * param.potential_scale - * Q_rxn - / param.L_x, - "X-averaged electrochemical heating": Q_rxn_av, - "X-averaged electrochemical heating [A.V.m-3]": param.i_typ - * param.potential_scale - * Q_rxn_av - / param.L_x, - "Volume-averaged electrochemical heating": Q_rxn_vol_av, - "Volume-averaged electrochemical heating [A.V.m-3]": param.i_typ - * param.potential_scale - * Q_rxn_vol_av - / param.L_x, + "Irreversible electrochemical heating [W.m-3]": Q_rxn * Q_scale, + "X-averaged irreversible electrochemical heating": Q_rxn_av, + "X-averaged irreversible electrochemical heating [W.m-3]": Q_rxn_av + * Q_scale, + "Volume-averaged irreversible electrochemical heating": Q_rxn_vol_av, + "Volume-averaged irreversible electrochemical heating [W.m-3]": Q_rxn_vol_av + * Q_scale, "Reversible heating": Q_rev, - "Reversible heating [A.V.m-3]": param.i_typ - * param.potential_scale - * Q_rev - / param.L_x, + "Reversible heating [W.m-3]": Q_rev * Q_scale, "X-averaged reversible heating": Q_rev_av, - "X-averaged reversible heating [A.V.m-3]": param.i_typ - * param.potential_scale - * Q_rev_av - / param.L_x, + "X-averaged reversible heating [W.m-3]": Q_rev_av * Q_scale, "Volume-averaged reversible heating": Q_rev_vol_av, - "Volume-averaged reversible heating [A.V.m-3]": param.i_typ - * param.potential_scale - * Q_rev_vol_av - / param.L_x, + "Volume-averaged reversible heating [W.m-3]": Q_rev_vol_av * Q_scale, "Total heating": Q, - "Total heating [A.V.m-3]": param.i_typ - * param.potential_scale - * Q - / param.L_x, + "Total heating [W.m-3]": Q * Q_scale, "X-averaged total heating": Q_av, - "X-averaged total heating [A.V.m-3]": param.i_typ - * param.potential_scale - * Q_av - / param.L_x, + "X-averaged total heating [W.m-3]": Q_av * Q_scale, "Volume-averaged total heating": Q_vol_av, - "Volume-averaged total heating [A.V.m-3]": param.i_typ - * param.potential_scale - * Q_vol_av - / param.L_x, + "Volume-averaged total heating [W.m-3]": Q_vol_av * Q_scale, } ) return variables + def _get_domain_heat_source_terms(self, Q_ohm_s, Q_ohm_e, Q_rxn, Q_rev): + """ + A private function to obtain the heat source terms split + by domain: 'negative electrode', 'separator' and 'positive electrode'. + + Parameters + ---------- + Q_ohm_s : :class:`pybamm.Symbol` + The Ohmic heating in the solid. + Q_ohm_e : :class:`pybamm.Symbol` + The Ohmic heating in the electrolyte. + Q_rxn : :class:`pybamm.Symbol` + The irreversible electrochemical heating. + Q_ohm_s : :class:`pybamm.Symbol` + The reversible electrochemical heating. + + Returns + ------- + variables : dict + The variables corresponding to the heat source terms in the + separate domains. + """ + param = self.param + + # Unpack by domain + Q_ohm_s_n, _, Q_ohm_s_p = Q_ohm_s.orphans + Q_ohm_e_n, Q_ohm_e_s, Q_ohm_e_p = Q_ohm_e.orphans + Q_ohm_n = Q_ohm_s_n + Q_ohm_e_n + Q_ohm_p = Q_ohm_s_p + Q_ohm_e_p + Q_rxn_n, _, Q_rxn_p = Q_rxn.orphans + Q_rev_n, _, Q_rev_p = Q_rev.orphans + Q_n = Q_ohm_n + Q_rxn_n + Q_rev_n + Q_s = Q_ohm_e_s + Q_p = Q_ohm_p + Q_rxn_p + Q_rev_p + + # X-average + Q_ohm_n_av = pybamm.x_average(Q_ohm_n) + Q_ohm_s_av = pybamm.x_average(Q_ohm_e_s) + Q_ohm_p_av = pybamm.x_average(Q_ohm_p) + Q_rxn_n_av = pybamm.x_average(Q_rxn_n) + Q_rxn_p_av = pybamm.x_average(Q_rxn_n) + Q_rev_n_av = pybamm.x_average(Q_rev_n) + Q_rev_p_av = pybamm.x_average(Q_rev_n) + Q_n_av = pybamm.x_average(Q_n) + Q_s_av = pybamm.x_average(Q_s) + Q_p_av = pybamm.x_average(Q_p) + + # Volume average + Q_ohm_n_vol_av = self._yz_average(Q_ohm_n_av) + Q_ohm_s_vol_av = self._yz_average(Q_ohm_s_av) + Q_ohm_p_vol_av = self._yz_average(Q_ohm_p_av) + Q_rxn_n_vol_av = self._yz_average(Q_rxn_n_av) + Q_rxn_p_vol_av = self._yz_average(Q_rxn_p_av) + Q_rev_n_vol_av = self._yz_average(Q_rev_n_av) + Q_rev_p_vol_av = self._yz_average(Q_rev_p_av) + Q_n_vol_av = self._yz_average(Q_n_av) + Q_s_vol_av = self._yz_average(Q_s_av) + Q_p_vol_av = self._yz_average(Q_p_av) + + # Dimensional scaling for heat source terms + Q_scale = param.i_typ * param.potential_scale / param.L_x + + variables = { + "Negative electrode Ohmic heating": Q_ohm_n, + "Negative electrode Ohmic heating [W.m-3]": Q_ohm_n * Q_scale, + "X-averaged negative electrode Ohmic heating": Q_ohm_n_av, + "X-averaged negative electrode Ohmic heating [W.m-3]": Q_ohm_n_av * Q_scale, + "Volume-averaged negative electrode ohm heating": Q_ohm_n_vol_av, + "Volume-averaged negative electrode ohm heating [W.m-3]": Q_ohm_n_vol_av + * Q_scale, + "Separator Ohmic heating": Q_ohm_s, + "Separator Ohmic heating [W.m-3]": Q_ohm_s * Q_scale, + "X-averaged separator Ohmic heating": Q_ohm_s_av, + "X-averaged separator Ohmic heating [W.m-3]": Q_ohm_s_av * Q_scale, + "Volume-averaged separator ohm heating": Q_ohm_s_vol_av, + "Volume-averaged separator ohm heating [W.m-3]": Q_ohm_s_vol_av * Q_scale, + "Positive electrode Ohmic heating": Q_ohm_p, + "Positive electrode Ohmic heating [W.m-3]": Q_ohm_p * Q_scale, + "X-averaged positive electrode Ohmic heating": Q_ohm_p_av, + "X-averaged positive electrode Ohmic heating [W.m-3]": Q_ohm_p_av * Q_scale, + "Volume-averaged positive electrode Ohmic heating": Q_ohm_p_vol_av, + "Volume-averaged positive electrode Ohmic heating [W.m-3]": Q_ohm_p_vol_av + * Q_scale, + "Negative electrode irreversible electrochemical heating": Q_rxn_n_av, + "Negative electrode irreversible electrochemical heating [W.m-3]": Q_rxn_n_av + * Q_scale, + "X-averaged negative electrode irreversible electrochemical heating": Q_rxn_n_av, + "X-averaged negative electrode irreversible electrochemical heating [W.m-3]": Q_rxn_n_av + * Q_scale, + "Volume-averaged negative electrode irreversible electrochemical heating": Q_rxn_n_vol_av, + "Volume-averaged negative electrode irreversible electrochemical heating [W.m-3]": Q_rxn_n_vol_av + * Q_scale, + "Positive electrode irreversible electrochemical heating": Q_rxn_p_av, + "Positive electrode irreversible electrochemical heating [W.m-3]": Q_rxn_p_av + * Q_scale, + "X-averaged positive electrode irreversible electrochemical heating": Q_rxn_p_av, + "X-averaged positive electrode irreversible electrochemical heating [W.m-3]": Q_rxn_p_av + * Q_scale, + "Volume-averaged positive electrode irreversible electrochemical heating": Q_rxn_p_vol_av, + "Volume-averaged positive electrode irreversible electrochemical heating [W.m-3]": Q_rxn_p_vol_av + * Q_scale, + "Negative electrode reversible heating": Q_rev_n, + "Negative electrode reversible heating [W.m-3]": Q_rev_n * Q_scale, + "X-averaged negative electrode reversible heating": Q_rev_n_av, + "X-averaged negative electrode reversible heating [W.m-3]": Q_rev_n_av + * Q_scale, + "Volume-averaged negative electrode reversible heating": Q_rev_n_vol_av, + "Volume-averaged negative electrode reversible heating [W.m-3]": Q_rev_n_vol_av + * Q_scale, + "Positive electrode reversible heating": Q_rev_p, + "Positive electrode reversible heating [W.m-3]": Q_rev_p * Q_scale, + "X-averaged positive electrode reversible heating": Q_rev_p_av, + "X-averaged positive electrode reversible heating [W.m-3]": Q_rev_p_av + * Q_scale, + "Volume-averaged positive electrode reversible heating": Q_rev_p_vol_av, + "Volume-averaged positive electrode reversible heating [W.m-3]": Q_rev_p_vol_av + * Q_scale, + "Negative electrode total heating": Q_n, + "Negative electrode total heating [W.m-3]": Q_n * Q_scale, + "X-averaged negative electrode total heating": Q_n_av, + "X-averaged negative electrode total heating [W.m-3]": Q_n_av * Q_scale, + "Volume-averaged negative electrode total heating": Q_n_vol_av, + "Volume-averaged negative electrode total heating [W.m-3]": Q_n_vol_av + * Q_scale, + "Separator total heating": Q_s, + "Separator total heating [W.m-3]": Q_s * Q_scale, + "X-averaged separator total heating": Q_s_av, + "X-averaged separator total heating [W.m-3]": Q_s_av * Q_scale, + "Volume-averaged separator total heating": Q_s_vol_av, + "Volume-averaged separator total heating [W.m-3]": Q_s_vol_av * Q_scale, + "Positive electrode total heating": Q_p, + "Positive electrode total heating [W.m-3]": Q_p * Q_scale, + "X-averaged postive electrode total heating": Q_p_av, + "X-averaged postive electrode total heating [W.m-3]": Q_p_av * Q_scale, + "Volume-averaged postive electrode total heating": Q_p_vol_av, + "Volume-averaged positive electrode total heating [W.m-3]": Q_p_vol_av + * Q_scale, + } + + return variables + def _flux_law(self, T): raise NotImplementedError @@ -252,3 +386,16 @@ def _x_average(self, var, var_cn, var_cp): + self.param.l_cp * var_cp ) / self.param.l return out + + def _effective_properties(self): + """ + Computes the effective effective product of density and specific heat, and + effective thermal conductivity, respectively. These are computed differently + depending upon whether current collectors are included or not. Defualt + behaviour is to assume the presence of current collectors. Due to the choice + of non-dimensionalisation, the dimensionless effective properties are equal + to 1 in the case where current collectors are accounted for. + """ + rho_eff = pybamm.Scalar(1) + lambda_eff = pybamm.Scalar(1) + return rho_eff, lambda_eff diff --git a/pybamm/models/submodels/thermal/isothermal/isothermal.py b/pybamm/models/submodels/thermal/isothermal/isothermal.py index 20dc2721ba..f8ecdfc5ab 100644 --- a/pybamm/models/submodels/thermal/isothermal/isothermal.py +++ b/pybamm/models/submodels/thermal/isothermal/isothermal.py @@ -39,17 +39,17 @@ def get_coupled_variables(self, variables): variables.update( { "Ohmic heating": pybamm.Scalar(0), - "Ohmic heating [A.V.m-3]": pybamm.Scalar(0), + "Ohmic heating [W.m-3]": pybamm.Scalar(0), "Irreversible electrochemical heating": pybamm.Scalar(0), - "Irreversible electrochemical heating [A.V.m-3]": pybamm.Scalar(0), + "Irreversible electrochemical heating [W.m-3]": pybamm.Scalar(0), "Reversible heating": pybamm.Scalar(0), - "Reversible heating [A.V.m-3]": pybamm.Scalar(0), + "Reversible heating [W.m-3]": pybamm.Scalar(0), "Total heating": pybamm.Scalar(0), - "Total heating [A.V.m-3]": pybamm.Scalar(0), + "Total heating [W.m-3]": pybamm.Scalar(0), "X-averaged total heating": pybamm.Scalar(0), - "X-averaged total heating [A.V.m-3]": pybamm.Scalar(0), + "X-averaged total heating [W.m-3]": pybamm.Scalar(0), "Volume-averaged total heating": pybamm.Scalar(0), - "Volume-averaged total heating [A.V.m-3]": pybamm.Scalar(0), + "Volume-averaged total heating [W.m-3]": pybamm.Scalar(0), } ) return variables diff --git a/pybamm/models/submodels/thermal/x_full/x_full_no_current_collector.py b/pybamm/models/submodels/thermal/x_full/x_full_no_current_collector.py index 9d81bf76ca..a919e6003b 100644 --- a/pybamm/models/submodels/thermal/x_full/x_full_no_current_collector.py +++ b/pybamm/models/submodels/thermal/x_full/x_full_no_current_collector.py @@ -50,8 +50,10 @@ def _current_collector_heating(self, variables): return Q_s_cn, Q_s_cp def _yz_average(self, var): - """Computes the y-z average by integration over y and z - In this case this is just equal to the input variable""" + """ + Computes the y-z average by integration over y and z + In this case this is just equal to the input variable + """ return var def _x_average(self, var, var_cn, var_cp): diff --git a/pybamm/models/submodels/thermal/x_lumped/base_x_lumped.py b/pybamm/models/submodels/thermal/x_lumped/base_x_lumped.py index 3114a2ceb2..fa93592e60 100644 --- a/pybamm/models/submodels/thermal/x_lumped/base_x_lumped.py +++ b/pybamm/models/submodels/thermal/x_lumped/base_x_lumped.py @@ -51,3 +51,11 @@ def _flux_law(self, T): def set_initial_conditions(self, variables): T = variables["X-averaged cell temperature"] self.initial_conditions = {T: self.param.T_init} + + def _surface_cooling_coefficient(self): + """Returns the surface cooling coefficient in for x-lumped models.""" + # Account for surface area to volume ratio in cooling coefficient + # Note: assumes pouch cell geometry + A = self.param.l_y * self.param.l_z + V = self.param.l * self.param.l_y * self.param.l_z + return -2 * self.param.h * A / V / (self.param.delta ** 2) diff --git a/pybamm/models/submodels/thermal/x_lumped/x_lumped_0D_current_collectors.py b/pybamm/models/submodels/thermal/x_lumped/x_lumped_0D_current_collectors.py index d705c0d684..1d2025a16a 100644 --- a/pybamm/models/submodels/thermal/x_lumped/x_lumped_0D_current_collectors.py +++ b/pybamm/models/submodels/thermal/x_lumped/x_lumped_0D_current_collectors.py @@ -14,12 +14,10 @@ def set_rhs(self, variables): T_av = variables["X-averaged cell temperature"] Q_av = variables["X-averaged total heating"] + cooling_coeff = self._surface_cooling_coefficient() + self.rhs = { - T_av: ( - self.param.B * Q_av - - (2 * self.param.h / (self.param.delta ** 2) / self.param.l) * T_av - ) - / self.param.C_th + T_av: (self.param.B * Q_av + cooling_coeff * T_av) / self.param.C_th } def _current_collector_heating(self, variables): diff --git a/pybamm/models/submodels/thermal/x_lumped/x_lumped_1D_current_collectors.py b/pybamm/models/submodels/thermal/x_lumped/x_lumped_1D_current_collectors.py index 71599b3203..fa96d98e70 100644 --- a/pybamm/models/submodels/thermal/x_lumped/x_lumped_1D_current_collectors.py +++ b/pybamm/models/submodels/thermal/x_lumped/x_lumped_1D_current_collectors.py @@ -15,12 +15,11 @@ def __init__(self, param): def set_rhs(self, variables): T_av = variables["X-averaged cell temperature"] Q_av = variables["X-averaged total heating"] + + cooling_coeff = self._surface_cooling_coefficient() + self.rhs = { - T_av: ( - pybamm.laplacian(T_av) - + self.param.B * Q_av - - (2 * self.param.h / (self.param.delta ** 2) / self.param.l) * T_av - ) + T_av: (pybamm.laplacian(T_av) + self.param.B * Q_av + cooling_coeff * T_av) / self.param.C_th } diff --git a/pybamm/models/submodels/thermal/x_lumped/x_lumped_2D_current_collectors.py b/pybamm/models/submodels/thermal/x_lumped/x_lumped_2D_current_collectors.py index a8117bd5b5..b9d301cc14 100644 --- a/pybamm/models/submodels/thermal/x_lumped/x_lumped_2D_current_collectors.py +++ b/pybamm/models/submodels/thermal/x_lumped/x_lumped_2D_current_collectors.py @@ -16,6 +16,8 @@ def set_rhs(self, variables): T_av = variables["X-averaged cell temperature"] Q_av = variables["X-averaged total heating"] + cooling_coeff = self._surface_cooling_coefficient() + # Add boundary source term which accounts for surface cooling around # the edge of the domain in the weak formulation. # TODO: update to allow different cooling conditions at the tabs @@ -23,8 +25,7 @@ def set_rhs(self, variables): T_av: ( pybamm.laplacian(T_av) + self.param.B * pybamm.source(Q_av, T_av) - - (2 * self.param.h / (self.param.delta ** 2) / self.param.l) - * pybamm.source(T_av, T_av) + + cooling_coeff * pybamm.source(T_av, T_av) - (self.param.h / self.param.delta) * pybamm.source(T_av, T_av, boundary=True) ) diff --git a/pybamm/models/submodels/thermal/x_lumped/x_lumped_no_current_collectors.py b/pybamm/models/submodels/thermal/x_lumped/x_lumped_no_current_collectors.py index 1e35fd4b88..db66a8f36f 100644 --- a/pybamm/models/submodels/thermal/x_lumped/x_lumped_no_current_collectors.py +++ b/pybamm/models/submodels/thermal/x_lumped/x_lumped_no_current_collectors.py @@ -28,12 +28,13 @@ def set_rhs(self, variables): T_av = variables["X-averaged cell temperature"] Q_av = variables["X-averaged total heating"] + # Get effective properties + rho_eff, _ = self._effective_properties() + cooling_coeff = self._surface_cooling_coefficient() + self.rhs = { - T_av: ( - self.param.B * Q_av - - (2 * self.param.h / (self.param.delta ** 2) / self.param.l) * T_av - ) - / self.param.C_th + T_av: (self.param.B * Q_av + cooling_coeff * T_av) + / (self.param.C_th * rho_eff) } def _current_collector_heating(self, variables): @@ -42,6 +43,18 @@ def _current_collector_heating(self, variables): Q_s_cp = pybamm.Scalar(0) return Q_s_cn, Q_s_cp + def _surface_cooling_coefficient(self): + """ + Returns the surface cooling coefficient in the absence of current + collectors. + """ + # Account for surface area to volume ratio in cooling coefficient + # Note: assumes pouch cell geometry and volume doesn't include current + # collectors + A = self.param.l_y * self.param.l_z + V = self.param.l_x * self.param.l_y * self.param.l_z + return -2 * self.param.h * A / V / (self.param.delta ** 2) + def _yz_average(self, var): """In 1D volume-averaged quantities are unchanged""" return var @@ -52,3 +65,26 @@ def _x_average(self, var, var_cn, var_cp): collectors. This overwrites the default behaviour of 'base_thermal'. """ return pybamm.x_average(var) + + def _effective_properties(self): + """ + Computes the effective effective product of density and specific heat, and + effective thermal conductivity, respectively. This overwrites the default + behaviour of 'base_thermal'. + """ + l_n = self.param.l_n + l_s = self.param.l_s + l_p = self.param.l_p + rho_n = self.param.rho_n + rho_s = self.param.rho_s + rho_p = self.param.rho_p + lambda_n = self.param.lambda_n + lambda_s = self.param.lambda_s + lambda_p = self.param.lambda_p + + rho_eff = (rho_n * l_n + rho_s * l_s + rho_p * l_p) / (l_n + l_n + l_p) + lambda_eff = (lambda_n * l_n + lambda_s * l_s + lambda_p * l_p) / ( + l_n + l_n + l_p + ) + + return rho_eff, lambda_eff diff --git a/pybamm/models/submodels/thermal/xyz_lumped/xyz_lumped_1D_current_collector.py b/pybamm/models/submodels/thermal/xyz_lumped/xyz_lumped_1D_current_collector.py index ffca162efa..a4919709f4 100644 --- a/pybamm/models/submodels/thermal/xyz_lumped/xyz_lumped_1D_current_collector.py +++ b/pybamm/models/submodels/thermal/xyz_lumped/xyz_lumped_1D_current_collector.py @@ -32,6 +32,7 @@ def _current_collector_heating(self, variables): def _surface_cooling_coefficient(self): """Returns the surface cooling coefficient in 1+1D""" + # Note: assumes pouch cell geometry return ( -2 * self.param.h / (self.param.delta ** 2) / self.param.l - self.param.l_z * self.param.h / self.param.delta diff --git a/pybamm/models/submodels/thermal/xyz_lumped/xyz_lumped_2D_current_collector.py b/pybamm/models/submodels/thermal/xyz_lumped/xyz_lumped_2D_current_collector.py index fcfbbd5329..bf2f9f3c5e 100644 --- a/pybamm/models/submodels/thermal/xyz_lumped/xyz_lumped_2D_current_collector.py +++ b/pybamm/models/submodels/thermal/xyz_lumped/xyz_lumped_2D_current_collector.py @@ -33,6 +33,7 @@ def _current_collector_heating(self, variables): def _surface_cooling_coefficient(self): """Returns the surface cooling coefficient in 2+1D""" + # Note: assumes pouch cell geometry return ( -2 * self.param.h / (self.param.delta ** 2) / self.param.l - 2 * (self.param.l_y + self.param.l_z) * self.param.h / self.param.delta diff --git a/pybamm/parameters/geometric_parameters.py b/pybamm/parameters/geometric_parameters.py index 1d7ea4780b..3b70db4865 100644 --- a/pybamm/parameters/geometric_parameters.py +++ b/pybamm/parameters/geometric_parameters.py @@ -53,11 +53,13 @@ l_s = L_s / L_x l_p = L_p / L_x l_cp = L_cp / L_x +l_x = L_x / L_x l_y = L_y / L_z l_z = L_z / L_z +a_cc = l_y * l_z l = L / L_x -delta = L_x / L_z +delta = L_x / L_z # Aspect ratio # Tab geometry l_tab_n = L_tab_n / L_z diff --git a/pybamm/parameters/standard_parameters_lithium_ion.py b/pybamm/parameters/standard_parameters_lithium_ion.py index 69403e3bd2..e981b597fd 100644 --- a/pybamm/parameters/standard_parameters_lithium_ion.py +++ b/pybamm/parameters/standard_parameters_lithium_ion.py @@ -38,7 +38,6 @@ L = pybamm.geometric_parameters.L A_cc = pybamm.geometric_parameters.A_cc - # Tab geometry L_tab_n = pybamm.geometric_parameters.L_tab_n Centre_y_tab_n = pybamm.geometric_parameters.Centre_y_tab_n @@ -254,8 +253,10 @@ def U_p_dimensional(sto, T): l_s = pybamm.geometric_parameters.l_s l_p = pybamm.geometric_parameters.l_p l_cp = pybamm.geometric_parameters.l_cp +l_x = pybamm.geometric_parameters.l_x l_y = pybamm.geometric_parameters.l_y l_z = pybamm.geometric_parameters.l_z +a_cc = pybamm.geometric_parameters.a_cc l = pybamm.geometric_parameters.l delta = pybamm.geometric_parameters.delta diff --git a/pybamm/parameters/thermal_parameters.py b/pybamm/parameters/thermal_parameters.py index 9c08124240..4041c4f5b1 100644 --- a/pybamm/parameters/thermal_parameters.py +++ b/pybamm/parameters/thermal_parameters.py @@ -38,12 +38,7 @@ "Positive current collector thermal conductivity [W.m-1.K-1]" ) -# Thermal parameters -h_dim = pybamm.Parameter("Heat transfer coefficient [W.m-2.K-1]") -Phi_dim = pybamm.Scalar(1) # typical scale for voltage drop across cell (order 1V) -Delta_T = ( - pybamm.electrical_parameters.i_typ * Phi_dim / h_dim -) # computed from balance of typical cross-cell Ohmic heating with surface heat loss +# Effective thermal properties rho_eff_dim = ( rho_cn_dim * c_p_cn_dim * pybamm.geometric_parameters.L_cn + rho_n_dim * c_p_n_dim * pybamm.geometric_parameters.L_n @@ -59,6 +54,15 @@ + lambda_cp_dim * pybamm.geometric_parameters.L_cp ) / pybamm.geometric_parameters.L +# Cooling coefficient +h_dim = pybamm.Parameter("Heat transfer coefficient [W.m-2.K-1]") + +# Typical temperature rise +Phi_dim = pybamm.Scalar(1) # typical scale for voltage drop across cell (order 1V) +Delta_T = ( + pybamm.electrical_parameters.i_typ * Phi_dim / h_dim +) # computed from balance of typical cross-cell Ohmic heating with surface heat loss + # Activation energies E_r_n = pybamm.Parameter("Negative reaction rate activation energy [J.mol-1]") E_r_p = pybamm.Parameter("Positive reaction rate activation energy [J.mol-1]") diff --git a/results/comsol_comparison/compare_comsol_thermal.py b/results/comsol_comparison/compare_comsol_thermal.py index 7c8d5d49a2..9c47f8053b 100644 --- a/results/comsol_comparison/compare_comsol_thermal.py +++ b/results/comsol_comparison/compare_comsol_thermal.py @@ -31,7 +31,7 @@ # create mesh var = pybamm.standard_spatial_vars -var_pts = {var.x_n: 31, var.x_s: 11, var.x_p: 31, var.r_n: 11, var.r_p: 11} +var_pts = {var.x_n: 101, var.x_s: 51, var.x_p: 101, var.r_n: 31, var.r_p: 31} mesh = pybamm.Mesh(geometry, pybamm_model.default_submesh_types, var_pts) # discretise model @@ -39,11 +39,14 @@ disc.process_model(pybamm_model) # discharge timescale -tau = param.process_symbol(pybamm.standard_parameters_lithium_ion.tau_discharge) +tau = param.process_symbol( + pybamm.standard_parameters_lithium_ion.tau_discharge +).evaluate(0) # solve model at comsol times -time = comsol_variables["time"] / tau.evaluate(0) -solution = pybamm_model.default_solver.solve(pybamm_model, time) +time = comsol_variables["time"] / tau +solver = pybamm.CasadiSolver() +solution = solver.solve(pybamm_model, time) "-----------------------------------------------------------------------------" "Make Comsol 'model' for comparison" @@ -88,6 +91,9 @@ def myinterp(t): comsol_phi_p = get_interp_fun(comsol_variables["phi_p"], ["positive electrode"]) comsol_voltage = interp.interp1d(comsol_t, comsol_variables["voltage"]) comsol_temperature = get_interp_fun(comsol_variables["temperature"], whole_cell) +comsol_temperature_av = interp.interp1d( + comsol_t, comsol_variables["average temperature"] +) comsol_q_irrev_n = get_interp_fun(comsol_variables["Q_irrev_n"], ["negative electrode"]) comsol_q_irrev_p = get_interp_fun(comsol_variables["Q_irrev_p"], ["positive electrode"]) comsol_q_rev_n = get_interp_fun(comsol_variables["Q_rev_n"], ["negative electrode"]) @@ -108,19 +114,202 @@ def myinterp(t): "Positive electrode potential [V]": comsol_phi_p, "Terminal voltage [V]": pybamm.Function(comsol_voltage, pybamm.t * tau), "Cell temperature [K]": comsol_temperature, + "Volume-averaged cell temperature [K]": pybamm.Function( + comsol_temperature_av, pybamm.t * tau + ), + "Negative electrode irreversible electrochemical heating [W.m-3]": comsol_q_irrev_n, + "Positive electrode irreversible electrochemical heating [W.m-3]": comsol_q_irrev_p, + "Negative electrode reversible heating [W.m-3]": comsol_q_rev_n, + "Positive electrode reversible heating [W.m-3]": comsol_q_rev_p, + "Negative electrode total heating [W.m-3]": comsol_q_total_n, + "Separator total heating [W.m-3]": comsol_q_total_s, + "Positive electrode total heating [W.m-3]": comsol_q_total_p, } "-----------------------------------------------------------------------------" "Plot comparison" -#TODO: fix QuickPlot -for var in comsol_model.variables.keys(): - plot = pybamm.QuickPlot( - [pybamm_model, comsol_model], - mesh, - [solution, solution], - output_variables=[var], - labels=["PyBaMM", "Comsol"], +plot_times = comsol_variables["time"] +pybamm_T = pybamm.ProcessedVariable( + pybamm_model.variables["Volume-averaged cell temperature [K]"], + solution.t, + solution.y, + mesh=mesh, +)(plot_times / tau) +comsol_T = pybamm.ProcessedVariable( + comsol_model.variables["Volume-averaged cell temperature [K]"], + solution.t, + solution.y, + mesh=mesh, +)(plot_times / tau) +plt.figure() +plt.plot(plot_times, pybamm_T, "-", label="PyBaMM") +plt.plot(plot_times, comsol_T, "o", label="COMSOL") +plt.xlabel("t") +plt.ylabel("T") +plt.legend() + +pybamm_voltage = pybamm.ProcessedVariable( + pybamm_model.variables["Terminal voltage [V]"], + solution.t, + solution.y, + mesh=mesh, +)(plot_times / tau) +comsol_voltage = pybamm.ProcessedVariable( + comsol_model.variables["Terminal voltage [V]"], + solution.t, + solution.y, + mesh=mesh, +)(plot_times / tau) +plt.figure() +plt.plot(plot_times, pybamm_voltage, "-", label="PyBaMM") +plt.plot(plot_times, comsol_voltage, "o", label="COMSOL") +plt.xlabel("t") +plt.ylabel("Voltage [V]") +plt.legend() + +# Get mesh nodes +x_n = mesh.combine_submeshes(*["negative electrode"])[0].nodes +x_s = mesh.combine_submeshes(*["separator"])[0].nodes +x_p = mesh.combine_submeshes(*["positive electrode"])[0].nodes +x = mesh.combine_submeshes(*whole_cell)[0].nodes + + +def comparison_plot(var, plot_times=None): + """ + Plot pybamm heat source (defined over whole cell) against comsol heat source + (defined by component) + """ + if plot_times is None: + plot_times = comsol_variables["time"] + + # Process pybamm heat source + pybamm_q = pybamm.ProcessedVariable( + pybamm_model.variables[var], solution.t, solution.y, mesh=mesh + ) + + # Process comsol heat source in negative electrode + comsol_q_n = pybamm.ProcessedVariable( + comsol_model.variables["Negative electrode " + var[0].lower() + var[1:]], + solution.t, + solution.y, + mesh=mesh, + ) + # Process comsol heat source in separator (if defined here) + try: + comsol_q_s = pybamm.ProcessedVariable( + comsol_model.variables["Separator " + var[0].lower() + var[1:]], + solution.t, + solution.y, + mesh=mesh, + ) + except KeyError: + comsol_q_s = None + print("Variable " + var + " not defined in separator") + # Process comsol heat source in positive electrode + comsol_q_p = pybamm.ProcessedVariable( + comsol_model.variables["Positive electrode " + var[0].lower() + var[1:]], + solution.t, + solution.y, + mesh=mesh, ) - plot.dynamic_plot(testing=True) + + # Make plot + if comsol_q_s: + n_cols = 3 + else: + n_cols = 2 + fig, ax = plt.subplots(1, n_cols, figsize=(15, 8)) + cmap = plt.get_cmap("inferno") + + for ind, t in enumerate(plot_times): + color = cmap(float(ind) / len(plot_times)) + ax[0].plot(x_n * L_x, pybamm_q(x=x_n, t=t / tau), "-", color=color) + ax[0].plot(x_n * L_x, comsol_q_n(x=x_n, t=t / tau), "o", color=color) + if comsol_q_s: + ax[1].plot(x_s * L_x, pybamm_q(x=x_s, t=t / tau), "-", color=color) + ax[1].plot(x_s * L_x, comsol_q_s(x=x_s, t=t / tau), "o", color=color) + ax[n_cols - 1].plot( + x_p * L_x, + pybamm_q(x=x_p, t=t / tau), + "-", + color=color, + label="PyBaMM" if ind == 0 else "", + ) + ax[n_cols - 1].plot( + x_p * L_x, + comsol_q_p(x=x_p, t=t / tau), + "o", + color=color, + label="COMSOL" if ind == 0 else "", + ) + + ax[0].set_xlabel("x_n") + ax[0].set_ylabel(var) + if comsol_q_s: + ax[1].set_xlabel("x_s") + ax[1].set_ylabel(var) + ax[n_cols - 1].set_xlabel("x_p") + ax[n_cols - 1].set_ylabel(var) + plt.legend() + plt.tight_layout() + + +def temperature_plot(plot_times=None): + """ + Plot pybamm heat source (defined over whole cell) against comsol heat source + (defined by component) + """ + if plot_times is None: + plot_times = comsol_variables["time"] + + # Process pybamm heat source + pybamm_T = pybamm.ProcessedVariable( + pybamm_model.variables["Cell temperature [K]"], + solution.t, + solution.y, + mesh=mesh, + ) + + # Process comsol heat source in negative electrode + comsol_T = pybamm.ProcessedVariable( + comsol_model.variables["Cell temperature [K]"], + solution.t, + solution.y, + mesh=mesh, + ) + + # Make plot + plt.figure(figsize=(15, 8)) + cmap = plt.get_cmap("inferno") + + for ind, t in enumerate(plot_times): + color = cmap(float(ind) / len(plot_times)) + plt.plot( + x * L_x, + pybamm_T(x=x, t=t / tau), + "-", + color=color, + label="PyBaMM" if ind == 0 else "", + ) + plt.plot( + x * L_x, + comsol_T(x=x, t=t / tau), + "o", + color=color, + label="COMSOL" if ind == 0 else "", + ) + + plt.xlabel("x") + plt.ylabel("Temperature [K]") + plt.legend() + plt.tight_layout() + + +# Make plots +plot_times = comsol_variables["time"][0::10] +comparison_plot("Irreversible electrochemical heating [W.m-3]", plot_times=plot_times) +comparison_plot("Reversible heating [W.m-3]", plot_times=plot_times) +comparison_plot("Total heating [W.m-3]", plot_times=plot_times) +# temperature_plot(plot_times) plt.show() diff --git a/results/comsol_comparison/load_comsol_data.py b/results/comsol_comparison/load_comsol_data.py index ea6f8366e0..a67ab9b8f3 100644 --- a/results/comsol_comparison/load_comsol_data.py +++ b/results/comsol_comparison/load_comsol_data.py @@ -65,6 +65,9 @@ comsol[1].values, (comsol_x_npts, comsol_time_npts), order="F" ) +# average temperature +comsol_temp_av = np.mean(comsol_temp_vals, axis=0) + # irreversible heat source in negative electrode comsol = pd.read_csv(path + "q_irrev_n.csv", sep=",", header=None) comsol_q_irrev_n_vals = np.reshape( @@ -125,6 +128,7 @@ "c_p_surf": comsol_c_p_surf_vals, "c_e": comsol_c_e_vals, "temperature": comsol_temp_vals, + "average temperature": comsol_temp_av, "Q_irrev_n": comsol_q_irrev_n_vals, "Q_irrev_p": comsol_q_irrev_p_vals, "Q_rev_n": comsol_q_rev_n_vals, From 2750870fa767d0cb5cfab5a744aeeabd10bdabb6 Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Wed, 30 Oct 2019 10:24:40 +0000 Subject: [PATCH 04/11] #678 adding missing electrolyte temp dependence --- .../base_electrolyte_conductivity.py | 2 +- ...igher_order_stefan_maxwell_conductivity.py | 33 +++++++++++-------- .../full_stefan_maxwell_conductivity.py | 3 +- ...urface_form_stefan_maxwell_conductivity.py | 16 ++++++--- 4 files changed, 35 insertions(+), 19 deletions(-) diff --git a/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py b/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py index cb3c49e64f..b282465c7b 100644 --- a/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py +++ b/pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py @@ -258,7 +258,7 @@ def _get_domain_current_variables(self, i_e, domain=None): def _get_whole_cell_variables(self, variables): """ A private function to obtain the potential and current concatenated - across the whole cell. Note required 'variables' to contain the potential + across the whole cell. Note: requires 'variables' to contain the potential and current in the subdomains: 'negative electrode', 'separator', and 'positive electrode'. diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_higher_order_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_higher_order_stefan_maxwell_conductivity.py index 0b482ba956..b6c5bf756c 100644 --- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_higher_order_stefan_maxwell_conductivity.py +++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_higher_order_stefan_maxwell_conductivity.py @@ -48,11 +48,10 @@ def get_coupled_variables(self, variables): eps_s_av = variables["Leading-order x-averaged separator porosity"] eps_p_av = variables["Leading-order x-averaged positive electrode porosity"] - # Note: here we want the average of the temperature over the negative - # electrode, separator and positive electrode (not including the current - # collectors) - T = variables["Cell temperature"] - T_av = pybamm.x_average(T) + T_av = variables["X-averaged cell temperature"] + T_av_n = pybamm.PrimaryBroadcast(T_av, "negative electrode") + T_av_s = pybamm.PrimaryBroadcast(T_av, "separator") + T_av_p = pybamm.PrimaryBroadcast(T_av, "positive electrode") c_e_n, c_e_s, c_e_p = c_e.orphans @@ -90,6 +89,7 @@ def get_coupled_variables(self, variables): + phi_s_n_av - ( chi_av + * (1 + param.Theta * T_av) * pybamm.x_average( self._higher_order_macinnes_function( c_e_n / pybamm.PrimaryBroadcast(c_e_av, "negative electrode") @@ -106,6 +106,7 @@ def get_coupled_variables(self, variables): pybamm.PrimaryBroadcast(phi_e_const, "negative electrode") + ( chi_av_n + * (1 + param.Theta * T_av_n) * self._higher_order_macinnes_function( c_e_n / pybamm.PrimaryBroadcast(c_e_av, "negative electrode") ) @@ -124,6 +125,7 @@ def get_coupled_variables(self, variables): pybamm.PrimaryBroadcast(phi_e_const, "separator") + ( chi_av_s + * (1 + param.Theta * T_av_s) * self._higher_order_macinnes_function( c_e_s / pybamm.PrimaryBroadcast(c_e_av, "separator") ) @@ -137,6 +139,7 @@ def get_coupled_variables(self, variables): pybamm.PrimaryBroadcast(phi_e_const, "positive electrode") + ( chi_av_p + * (1 + param.Theta * T_av_p) * self._higher_order_macinnes_function( c_e_p / pybamm.PrimaryBroadcast(c_e_av, "positive electrode") ) @@ -155,15 +158,19 @@ def get_coupled_variables(self, variables): phi_e_av = pybamm.x_average(phi_e) # concentration overpotential - eta_c_av = chi_av * ( - pybamm.x_average( - self._higher_order_macinnes_function( - c_e_p / pybamm.PrimaryBroadcast(c_e_av, "positive electrode") + eta_c_av = ( + chi_av + * (1 + param.Theta * T_av) + * ( + pybamm.x_average( + self._higher_order_macinnes_function( + c_e_p / pybamm.PrimaryBroadcast(c_e_av, "positive electrode") + ) ) - ) - - pybamm.x_average( - self._higher_order_macinnes_function( - c_e_n / pybamm.PrimaryBroadcast(c_e_av, "negative electrode") + - pybamm.x_average( + self._higher_order_macinnes_function( + c_e_n / pybamm.PrimaryBroadcast(c_e_av, "negative electrode") + ) ) ) ) diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py index 536c989094..f4fae2f790 100644 --- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py +++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py @@ -39,7 +39,8 @@ def get_coupled_variables(self, variables): phi_e = variables["Electrolyte potential"] i_e = (param.kappa_e(c_e, T) * (eps ** param.b) * param.gamma_e / param.C_e) * ( - param.chi(c_e) * pybamm.grad(c_e) / c_e - pybamm.grad(phi_e) + param.chi(c_e) * (1 + param.Theta * T) * pybamm.grad(c_e) / c_e + - pybamm.grad(phi_e) ) variables.update(self._get_standard_current_variables(i_e)) diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py index 5059986287..0c49892b9a 100644 --- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py +++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py @@ -73,13 +73,17 @@ def set_boundary_conditions(self, variables): i_boundary_cc = variables["Current collector current density"] c_e = variables[self.domain + " electrolyte concentration"] delta_phi = variables[self.domain + " electrode surface potential difference"] + T = variables["Cell temperature"] if self.domain == "Negative": c_e_flux = pybamm.BoundaryGradient(c_e, "right") flux_left = -i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "left") flux_right = ( (i_boundary_cc / pybamm.BoundaryValue(conductivity, "right")) - - pybamm.BoundaryValue(param.chi(c_e) / c_e, "right") * c_e_flux + - pybamm.BoundaryValue( + (1 + param.Theta * T) * param.chi(c_e) / c_e, "right" + ) + * c_e_flux - i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "right") ) @@ -92,7 +96,10 @@ def set_boundary_conditions(self, variables): c_e_flux = pybamm.BoundaryGradient(c_e, "left") flux_left = ( (i_boundary_cc / pybamm.BoundaryValue(conductivity, "left")) - - pybamm.BoundaryValue(param.chi(c_e) / c_e, "left") * c_e_flux + - pybamm.BoundaryValue( + (1 + param.Theta * T) * param.chi(c_e) / c_e, "left" + ) + * c_e_flux - i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "left") ) flux_right = -i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "right") @@ -152,9 +159,10 @@ def _get_neg_pos_coupled_variables(self, variables): i_boundary_cc = variables["Current collector current density"] c_e = variables[self.domain + " electrolyte concentration"] delta_phi = variables[self.domain + " electrode surface potential difference"] + T = variables[self.domain + " temperature"] i_e = conductivity * ( - (param.chi(c_e) / c_e) * pybamm.grad(c_e) + ((1 + param.Theta * T) * param.chi(c_e) / c_e) * pybamm.grad(c_e) + pybamm.grad(delta_phi) + pybamm.PrimaryBroadcast(i_boundary_cc, self.domain_for_broadcast) / sigma_eff @@ -190,7 +198,7 @@ def _get_sep_coupled_variables(self, variables): phi_e_s = pybamm.PrimaryBroadcast( pybamm.boundary_value(phi_e_n, "right"), "separator" ) + pybamm.IndefiniteIntegral( - chi_e_s / c_e_s * pybamm.grad(c_e_s) + (1 + param.Theta * T) * chi_e_s / c_e_s * pybamm.grad(c_e_s) - param.C_e * pybamm.PrimaryBroadcast(i_boundary_cc, self.domain_for_broadcast) / kappa_s_eff, From 57c73c7a09ae976a273d6222a12184d0f0084e2b Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Wed, 30 Oct 2019 11:11:19 +0000 Subject: [PATCH 05/11] #678 debug electrolyte --- .../full_surface_form_stefan_maxwell_conductivity.py | 5 +++-- pybamm/models/submodels/thermal/base_thermal.py | 4 ++-- pybamm/models/submodels/thermal/isothermal/isothermal.py | 4 ++-- .../test_composite_stefan_maxwell_conductivity.py | 2 +- .../test_full_surface_form_stefan_maxwell_conductivity.py | 2 ++ 5 files changed, 10 insertions(+), 7 deletions(-) diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py index 0c49892b9a..996f9fdfe3 100644 --- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py +++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py @@ -73,9 +73,9 @@ def set_boundary_conditions(self, variables): i_boundary_cc = variables["Current collector current density"] c_e = variables[self.domain + " electrolyte concentration"] delta_phi = variables[self.domain + " electrode surface potential difference"] - T = variables["Cell temperature"] if self.domain == "Negative": + T = variables["Negative electrode temperature"] c_e_flux = pybamm.BoundaryGradient(c_e, "right") flux_left = -i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "left") flux_right = ( @@ -93,6 +93,7 @@ def set_boundary_conditions(self, variables): rbc_c_e = (c_e_flux, "Neumann") elif self.domain == "Positive": + T = variables["Positive electrode temperature"] c_e_flux = pybamm.BoundaryGradient(c_e, "left") flux_left = ( (i_boundary_cc / pybamm.BoundaryValue(conductivity, "left")) @@ -159,7 +160,7 @@ def _get_neg_pos_coupled_variables(self, variables): i_boundary_cc = variables["Current collector current density"] c_e = variables[self.domain + " electrolyte concentration"] delta_phi = variables[self.domain + " electrode surface potential difference"] - T = variables[self.domain + " temperature"] + T = variables[self.domain + " electrode temperature"] i_e = conductivity * ( ((1 + param.Theta * T) * param.chi(c_e) / c_e) * pybamm.grad(c_e) diff --git a/pybamm/models/submodels/thermal/base_thermal.py b/pybamm/models/submodels/thermal/base_thermal.py index 4495694c82..f2b16de228 100644 --- a/pybamm/models/submodels/thermal/base_thermal.py +++ b/pybamm/models/submodels/thermal/base_thermal.py @@ -179,8 +179,8 @@ def _get_standard_coupled_variables(self, variables): "X-averaged irreversible electrochemical heating [W.m-3]": Q_rxn_av * Q_scale, "Volume-averaged irreversible electrochemical heating": Q_rxn_vol_av, - "Volume-averaged irreversible electrochemical heating [W.m-3]": Q_rxn_vol_av - * Q_scale, + "Volume-averaged irreversible electrochemical heating [W.m-3]": + Q_rxn_vol_av * Q_scale, "Reversible heating": Q_rev, "Reversible heating [W.m-3]": Q_rev * Q_scale, "X-averaged reversible heating": Q_rev_av, diff --git a/pybamm/models/submodels/thermal/isothermal/isothermal.py b/pybamm/models/submodels/thermal/isothermal/isothermal.py index f8ecdfc5ab..3c12d58d5e 100644 --- a/pybamm/models/submodels/thermal/isothermal/isothermal.py +++ b/pybamm/models/submodels/thermal/isothermal/isothermal.py @@ -73,6 +73,6 @@ def _yz_average(self, var): def _x_average(self, var, var_cn, var_cp): """ Temperature is uniform and heat source terms are zero, so the average - returns the input variable. + returns zeros broadcasted onto the current collector domain. """ - return var + return pybamm.PrimaryBroadcast(0, "current collector") diff --git a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_composite_stefan_maxwell_conductivity.py b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_composite_stefan_maxwell_conductivity.py index 6eff2e9d01..c3e55d659f 100644 --- a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_composite_stefan_maxwell_conductivity.py +++ b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_composite_stefan_maxwell_conductivity.py @@ -22,7 +22,7 @@ def test_public_functions(self): "Leading-order x-averaged negative electrode porosity": a, "Leading-order x-averaged separator porosity": a, "Leading-order x-averaged positive electrode porosity": a, - "Cell temperature": a, + "X-averaged cell temperature": a, } submodel = pybamm.electrolyte.stefan_maxwell.conductivity.Composite(param) std_tests = tests.StandardSubModelTests(submodel, variables) diff --git a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py index 464cb91915..bc4bc47d73 100644 --- a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py +++ b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py @@ -27,6 +27,8 @@ def test_public_functions(self): "Negative electrode interfacial current density": a_n, "Electrolyte potential": pybamm.Concatenation(a_n, a_s, a_p), "Negative electrode temperature": a_n, + "Separator temperature": a_s, + "Positive electrode temperature": a_p, } icd = " interfacial current density" reactions = { From 9940bf76a6a3a0daad95ef114243c02b575fe71d Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Wed, 30 Oct 2019 14:45:10 +0000 Subject: [PATCH 06/11] #678 fix temp depend kinetics --- .../inverse_kinetics/base_inverse_kinetics.py | 9 +- .../inverse_kinetics/inverse_butler_volmer.py | 6 +- .../interface/kinetics/base_kinetics.py | 25 ++- .../interface/kinetics/butler_volmer.py | 19 ++- .../interface/kinetics/no_reaction.py | 2 +- .../submodels/interface/kinetics/tafel.py | 23 ++- .../models/submodels/thermal/base_thermal.py | 147 ------------------ .../xyz_lumped_2D_current_collector.py | 2 +- .../test_interface/test_lead_acid.py | 2 + 9 files changed, 59 insertions(+), 176 deletions(-) diff --git a/pybamm/models/submodels/interface/inverse_kinetics/base_inverse_kinetics.py b/pybamm/models/submodels/interface/inverse_kinetics/base_inverse_kinetics.py index 6c7acfe61c..244f3d2d65 100644 --- a/pybamm/models/submodels/interface/inverse_kinetics/base_inverse_kinetics.py +++ b/pybamm/models/submodels/interface/inverse_kinetics/base_inverse_kinetics.py @@ -41,8 +41,13 @@ def get_coupled_variables(self, variables): ne = self.param.ne_n elif self.domain == "Positive": ne = self.param.ne_p + # Note: T must have the same domain as j0 and eta_r + if j0.domain in ["current collector", ["current collector"]]: + T = variables["X-averaged cell temperature"] + else: + T = variables[self.domain + " electrode temperature"] - eta_r = self._get_overpotential(j, j0, ne) + eta_r = self._get_overpotential(j, j0, ne, T) delta_phi = eta_r + ocp variables.update(self._get_standard_interfacial_current_variables(j)) @@ -73,5 +78,5 @@ def get_coupled_variables(self, variables): return variables - def _get_overpotential(self, j, j0, ne): + def _get_overpotential(self, j, j0, ne, T): raise NotImplementedError diff --git a/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py b/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py index 37afd228ab..7edd3b5753 100644 --- a/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py +++ b/pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py @@ -28,5 +28,7 @@ class InverseButlerVolmer(BaseInverseKinetics, ButlerVolmer): def __init__(self, param, domain): super().__init__(param, domain) - def _get_overpotential(self, j, j0, ne): - return (2 / ne) * pybamm.Function(np.arcsinh, j / (2 * j0)) + def _get_overpotential(self, j, j0, ne, T): + return (2 * (1 + self.param.Theta * T) / ne) * pybamm.Function( + np.arcsinh, j / (2 * j0) + ) diff --git a/pybamm/models/submodels/interface/kinetics/base_kinetics.py b/pybamm/models/submodels/interface/kinetics/base_kinetics.py index a4cfe9f33e..b65181e792 100644 --- a/pybamm/models/submodels/interface/kinetics/base_kinetics.py +++ b/pybamm/models/submodels/interface/kinetics/base_kinetics.py @@ -40,8 +40,13 @@ def get_coupled_variables(self, variables): eta_r = delta_phi - ocp # Get number of electrons in reaction ne = self._get_number_of_electrons_in_reaction() - - j = self._get_kinetics(j0, ne, eta_r) + # Get kinetics. Note: T must have the same domain as j0 and eta_r + if j0.domain in ["current collector", ["current collector"]]: + T = variables["X-averaged cell temperature"] + else: + T = variables[self.domain + " electrode temperature"] + j = self._get_kinetics(j0, ne, eta_r, T) + # Get average interfacial current density j_tot_av = self._get_average_total_interfacial_current_density(variables) # j = j_tot_av + (j - pybamm.x_average(j)) # enforce true average @@ -73,7 +78,7 @@ def get_coupled_variables(self, variables): def _get_exchange_current_density(self, variables): raise NotImplementedError - def _get_kinetics(self, j0, ne, eta_r): + def _get_kinetics(self, j0, ne, eta_r, T): raise NotImplementedError def _get_open_circuit_potential(self, variables): @@ -84,10 +89,10 @@ def _get_dj_dc(self, variables): Default to calculate derivative of interfacial current density with respect to concentration. Can be overwritten by specific kinetic functions. """ - c_e, delta_phi, j0, ne, ocp = self._get_interface_variables_for_first_order( + c_e, delta_phi, j0, ne, ocp, T = self._get_interface_variables_for_first_order( variables ) - j = self._get_kinetics(j0, ne, delta_phi - ocp) + j = self._get_kinetics(j0, ne, delta_phi - ocp, T) return j.diff(c_e) def _get_dj_ddeltaphi(self, variables): @@ -95,10 +100,10 @@ def _get_dj_ddeltaphi(self, variables): Default to calculate derivative of interfacial current density with respect to surface potential difference. Can be overwritten by specific kinetic functions. """ - _, delta_phi, j0, ne, ocp = self._get_interface_variables_for_first_order( + _, delta_phi, j0, ne, ocp, T = self._get_interface_variables_for_first_order( variables ) - j = self._get_kinetics(j0, ne, delta_phi - ocp) + j = self._get_kinetics(j0, ne, delta_phi - ocp, T) return j.diff(delta_phi) def _get_interface_variables_for_first_order(self, variables): @@ -118,7 +123,11 @@ def _get_interface_variables_for_first_order(self, variables): j0 = self._get_exchange_current_density(hacked_variables) ne = self._get_number_of_electrons_in_reaction() ocp = self._get_open_circuit_potential(hacked_variables)[0] - return c_e_0, delta_phi, j0, ne, ocp + if j0.domain in ["current collector", ["current collector"]]: + T = variables["X-averaged cell temperature"] + else: + T = variables[self.domain + " electrode temperature"] + return c_e_0, delta_phi, j0, ne, ocp, T def _get_j_diffusion_limited_first_order(self, variables): """ diff --git a/pybamm/models/submodels/interface/kinetics/butler_volmer.py b/pybamm/models/submodels/interface/kinetics/butler_volmer.py index 437298bb84..0ee18a8f6b 100644 --- a/pybamm/models/submodels/interface/kinetics/butler_volmer.py +++ b/pybamm/models/submodels/interface/kinetics/butler_volmer.py @@ -12,7 +12,7 @@ class ButlerVolmer(BaseModel): Base submodel which implements the forward Butler-Volmer equation: .. math:: - j = j_0(c) * \\sinh(\\eta_r(c)) + j = 2 * j_0(c) * \\sinh( (ne / (2 * (1 + \\Theta T)) * \\eta_r(c)) Parameters ---------- @@ -28,26 +28,29 @@ class ButlerVolmer(BaseModel): def __init__(self, param, domain): super().__init__(param, domain) - def _get_kinetics(self, j0, ne, eta_r): - return 2 * j0 * pybamm.sinh((ne / 2) * eta_r) + def _get_kinetics(self, j0, ne, eta_r, T): + prefactor = ne / (2 * (1 + self.param.Theta * T)) + return 2 * j0 * pybamm.sinh(prefactor * eta_r) def _get_dj_dc(self, variables): "See :meth:`pybamm.interface.kinetics.BaseModel._get_dj_dc`" - c_e, delta_phi, j0, ne, ocp = self._get_interface_variables_for_first_order( + c_e, delta_phi, j0, ne, ocp, T = self._get_interface_variables_for_first_order( variables ) eta_r = delta_phi - ocp - return (2 * j0.diff(c_e) * pybamm.sinh((ne / 2) * eta_r)) - ( - 2 * j0 * (ne / 2) * ocp.diff(c_e) * pybamm.cosh((ne / 2) * eta_r) + prefactor = ne / (2 * (1 + self.param.Theta * T)) + return (2 * j0.diff(c_e) * pybamm.sinh(prefactor * eta_r)) - ( + 2 * j0 * prefactor * ocp.diff(c_e) * pybamm.cosh(prefactor * eta_r) ) def _get_dj_ddeltaphi(self, variables): "See :meth:`pybamm.interface.kinetics.BaseModel._get_dj_ddeltaphi`" - _, delta_phi, j0, ne, ocp = self._get_interface_variables_for_first_order( + _, delta_phi, j0, ne, ocp, T = self._get_interface_variables_for_first_order( variables ) eta_r = delta_phi - ocp - return 2 * j0 * (ne / 2) * pybamm.cosh((ne / 2) * eta_r) + prefactor = ne / (2 * (1 + self.param.Theta * T)) + return 2 * j0 * prefactor * pybamm.cosh(prefactor * eta_r) class FirstOrderButlerVolmer(ButlerVolmer, BaseFirstOrderKinetics): diff --git a/pybamm/models/submodels/interface/kinetics/no_reaction.py b/pybamm/models/submodels/interface/kinetics/no_reaction.py index e135c5f25c..d16c32a933 100644 --- a/pybamm/models/submodels/interface/kinetics/no_reaction.py +++ b/pybamm/models/submodels/interface/kinetics/no_reaction.py @@ -24,5 +24,5 @@ class NoReaction(BaseModel): def __init__(self, param, domain): super().__init__(param, domain) - def _get_kinetics(self, j0, ne, eta_r): + def _get_kinetics(self, j0, ne, eta_r, T): return pybamm.Scalar(0) diff --git a/pybamm/models/submodels/interface/kinetics/tafel.py b/pybamm/models/submodels/interface/kinetics/tafel.py index fbd0c5f770..a263657353 100644 --- a/pybamm/models/submodels/interface/kinetics/tafel.py +++ b/pybamm/models/submodels/interface/kinetics/tafel.py @@ -11,7 +11,7 @@ class ForwardTafel(BaseModel): Base submodel which implements the forward Tafel equation: .. math:: - j = j_0(c) * \\exp(\\eta_r(c)) + j = j_0(c) * \\exp((ne / (2 * (1 + \\Theta T)) * \\eta_r(c)) Parameters ---------- @@ -27,24 +27,33 @@ class ForwardTafel(BaseModel): def __init__(self, param, domain): super().__init__(param, domain) - def _get_kinetics(self, j0, ne, eta_r): - return j0 * pybamm.exp((ne / 2) * eta_r) + def _get_kinetics(self, j0, ne, eta_r, T): + return j0 * pybamm.exp((ne / (2 * (1 + self.param.Theta * T))) * eta_r) def _get_dj_dc(self, variables): "See :meth:`pybamm.interface.kinetics.BaseKinetics._get_dj_dc`" - c_e, delta_phi, j0, ne, ocp = self._get_interface_variables_for_first_order( + c_e, delta_phi, j0, ne, ocp, T = self._get_interface_variables_for_first_order( variables ) eta_r = delta_phi - ocp - return 2 * j0.diff(c_e) * pybamm.exp((ne / 2) * eta_r) + return ( + 2 + * j0.diff(c_e) + * pybamm.exp((ne / (2 * (1 + self.param.Theta * T))) * eta_r) + ) def _get_dj_ddeltaphi(self, variables): "See :meth:`pybamm.interface.kinetics.BaseKinetics._get_dj_ddeltaphi`" - _, delta_phi, j0, ne, ocp = self._get_interface_variables_for_first_order( + _, delta_phi, j0, ne, ocp, T = self._get_interface_variables_for_first_order( variables ) eta_r = delta_phi - ocp - return 2 * j0 * (ne / 2) * pybamm.exp((ne / 2) * eta_r) + return ( + 2 + * j0 + * (ne / (2 * (1 + self.param.Theta * T))) + * pybamm.exp((ne / 2) * eta_r) + ) class FirstOrderForwardTafel(ForwardTafel, BaseFirstOrderKinetics): diff --git a/pybamm/models/submodels/thermal/base_thermal.py b/pybamm/models/submodels/thermal/base_thermal.py index f2b16de228..2f7186b187 100644 --- a/pybamm/models/submodels/thermal/base_thermal.py +++ b/pybamm/models/submodels/thermal/base_thermal.py @@ -159,12 +159,6 @@ def _get_standard_coupled_variables(self, variables): # Dimensional scaling for heat source terms Q_scale = param.i_typ * param.potential_scale / param.L_x - # Geat heat source by domain - if isinstance(Q_ohm_e, pybamm.Concatenation): - variables.update( - self._get_domain_heat_source_terms(Q_ohm_s, Q_ohm_e, Q_rxn, Q_rev) - ) - variables.update( { "Ohmic heating": Q_ohm, @@ -197,147 +191,6 @@ def _get_standard_coupled_variables(self, variables): ) return variables - def _get_domain_heat_source_terms(self, Q_ohm_s, Q_ohm_e, Q_rxn, Q_rev): - """ - A private function to obtain the heat source terms split - by domain: 'negative electrode', 'separator' and 'positive electrode'. - - Parameters - ---------- - Q_ohm_s : :class:`pybamm.Symbol` - The Ohmic heating in the solid. - Q_ohm_e : :class:`pybamm.Symbol` - The Ohmic heating in the electrolyte. - Q_rxn : :class:`pybamm.Symbol` - The irreversible electrochemical heating. - Q_ohm_s : :class:`pybamm.Symbol` - The reversible electrochemical heating. - - Returns - ------- - variables : dict - The variables corresponding to the heat source terms in the - separate domains. - """ - param = self.param - - # Unpack by domain - Q_ohm_s_n, _, Q_ohm_s_p = Q_ohm_s.orphans - Q_ohm_e_n, Q_ohm_e_s, Q_ohm_e_p = Q_ohm_e.orphans - Q_ohm_n = Q_ohm_s_n + Q_ohm_e_n - Q_ohm_p = Q_ohm_s_p + Q_ohm_e_p - Q_rxn_n, _, Q_rxn_p = Q_rxn.orphans - Q_rev_n, _, Q_rev_p = Q_rev.orphans - Q_n = Q_ohm_n + Q_rxn_n + Q_rev_n - Q_s = Q_ohm_e_s - Q_p = Q_ohm_p + Q_rxn_p + Q_rev_p - - # X-average - Q_ohm_n_av = pybamm.x_average(Q_ohm_n) - Q_ohm_s_av = pybamm.x_average(Q_ohm_e_s) - Q_ohm_p_av = pybamm.x_average(Q_ohm_p) - Q_rxn_n_av = pybamm.x_average(Q_rxn_n) - Q_rxn_p_av = pybamm.x_average(Q_rxn_n) - Q_rev_n_av = pybamm.x_average(Q_rev_n) - Q_rev_p_av = pybamm.x_average(Q_rev_n) - Q_n_av = pybamm.x_average(Q_n) - Q_s_av = pybamm.x_average(Q_s) - Q_p_av = pybamm.x_average(Q_p) - - # Volume average - Q_ohm_n_vol_av = self._yz_average(Q_ohm_n_av) - Q_ohm_s_vol_av = self._yz_average(Q_ohm_s_av) - Q_ohm_p_vol_av = self._yz_average(Q_ohm_p_av) - Q_rxn_n_vol_av = self._yz_average(Q_rxn_n_av) - Q_rxn_p_vol_av = self._yz_average(Q_rxn_p_av) - Q_rev_n_vol_av = self._yz_average(Q_rev_n_av) - Q_rev_p_vol_av = self._yz_average(Q_rev_p_av) - Q_n_vol_av = self._yz_average(Q_n_av) - Q_s_vol_av = self._yz_average(Q_s_av) - Q_p_vol_av = self._yz_average(Q_p_av) - - # Dimensional scaling for heat source terms - Q_scale = param.i_typ * param.potential_scale / param.L_x - - variables = { - "Negative electrode Ohmic heating": Q_ohm_n, - "Negative electrode Ohmic heating [W.m-3]": Q_ohm_n * Q_scale, - "X-averaged negative electrode Ohmic heating": Q_ohm_n_av, - "X-averaged negative electrode Ohmic heating [W.m-3]": Q_ohm_n_av * Q_scale, - "Volume-averaged negative electrode ohm heating": Q_ohm_n_vol_av, - "Volume-averaged negative electrode ohm heating [W.m-3]": Q_ohm_n_vol_av - * Q_scale, - "Separator Ohmic heating": Q_ohm_s, - "Separator Ohmic heating [W.m-3]": Q_ohm_s * Q_scale, - "X-averaged separator Ohmic heating": Q_ohm_s_av, - "X-averaged separator Ohmic heating [W.m-3]": Q_ohm_s_av * Q_scale, - "Volume-averaged separator ohm heating": Q_ohm_s_vol_av, - "Volume-averaged separator ohm heating [W.m-3]": Q_ohm_s_vol_av * Q_scale, - "Positive electrode Ohmic heating": Q_ohm_p, - "Positive electrode Ohmic heating [W.m-3]": Q_ohm_p * Q_scale, - "X-averaged positive electrode Ohmic heating": Q_ohm_p_av, - "X-averaged positive electrode Ohmic heating [W.m-3]": Q_ohm_p_av * Q_scale, - "Volume-averaged positive electrode Ohmic heating": Q_ohm_p_vol_av, - "Volume-averaged positive electrode Ohmic heating [W.m-3]": Q_ohm_p_vol_av - * Q_scale, - "Negative electrode irreversible electrochemical heating": Q_rxn_n_av, - "Negative electrode irreversible electrochemical heating [W.m-3]": Q_rxn_n_av - * Q_scale, - "X-averaged negative electrode irreversible electrochemical heating": Q_rxn_n_av, - "X-averaged negative electrode irreversible electrochemical heating [W.m-3]": Q_rxn_n_av - * Q_scale, - "Volume-averaged negative electrode irreversible electrochemical heating": Q_rxn_n_vol_av, - "Volume-averaged negative electrode irreversible electrochemical heating [W.m-3]": Q_rxn_n_vol_av - * Q_scale, - "Positive electrode irreversible electrochemical heating": Q_rxn_p_av, - "Positive electrode irreversible electrochemical heating [W.m-3]": Q_rxn_p_av - * Q_scale, - "X-averaged positive electrode irreversible electrochemical heating": Q_rxn_p_av, - "X-averaged positive electrode irreversible electrochemical heating [W.m-3]": Q_rxn_p_av - * Q_scale, - "Volume-averaged positive electrode irreversible electrochemical heating": Q_rxn_p_vol_av, - "Volume-averaged positive electrode irreversible electrochemical heating [W.m-3]": Q_rxn_p_vol_av - * Q_scale, - "Negative electrode reversible heating": Q_rev_n, - "Negative electrode reversible heating [W.m-3]": Q_rev_n * Q_scale, - "X-averaged negative electrode reversible heating": Q_rev_n_av, - "X-averaged negative electrode reversible heating [W.m-3]": Q_rev_n_av - * Q_scale, - "Volume-averaged negative electrode reversible heating": Q_rev_n_vol_av, - "Volume-averaged negative electrode reversible heating [W.m-3]": Q_rev_n_vol_av - * Q_scale, - "Positive electrode reversible heating": Q_rev_p, - "Positive electrode reversible heating [W.m-3]": Q_rev_p * Q_scale, - "X-averaged positive electrode reversible heating": Q_rev_p_av, - "X-averaged positive electrode reversible heating [W.m-3]": Q_rev_p_av - * Q_scale, - "Volume-averaged positive electrode reversible heating": Q_rev_p_vol_av, - "Volume-averaged positive electrode reversible heating [W.m-3]": Q_rev_p_vol_av - * Q_scale, - "Negative electrode total heating": Q_n, - "Negative electrode total heating [W.m-3]": Q_n * Q_scale, - "X-averaged negative electrode total heating": Q_n_av, - "X-averaged negative electrode total heating [W.m-3]": Q_n_av * Q_scale, - "Volume-averaged negative electrode total heating": Q_n_vol_av, - "Volume-averaged negative electrode total heating [W.m-3]": Q_n_vol_av - * Q_scale, - "Separator total heating": Q_s, - "Separator total heating [W.m-3]": Q_s * Q_scale, - "X-averaged separator total heating": Q_s_av, - "X-averaged separator total heating [W.m-3]": Q_s_av * Q_scale, - "Volume-averaged separator total heating": Q_s_vol_av, - "Volume-averaged separator total heating [W.m-3]": Q_s_vol_av * Q_scale, - "Positive electrode total heating": Q_p, - "Positive electrode total heating [W.m-3]": Q_p * Q_scale, - "X-averaged postive electrode total heating": Q_p_av, - "X-averaged postive electrode total heating [W.m-3]": Q_p_av * Q_scale, - "Volume-averaged postive electrode total heating": Q_p_vol_av, - "Volume-averaged positive electrode total heating [W.m-3]": Q_p_vol_av - * Q_scale, - } - - return variables - def _flux_law(self, T): raise NotImplementedError diff --git a/pybamm/models/submodels/thermal/xyz_lumped/xyz_lumped_2D_current_collector.py b/pybamm/models/submodels/thermal/xyz_lumped/xyz_lumped_2D_current_collector.py index bf2f9f3c5e..60fe551a74 100644 --- a/pybamm/models/submodels/thermal/xyz_lumped/xyz_lumped_2D_current_collector.py +++ b/pybamm/models/submodels/thermal/xyz_lumped/xyz_lumped_2D_current_collector.py @@ -33,7 +33,7 @@ def _current_collector_heating(self, variables): def _surface_cooling_coefficient(self): """Returns the surface cooling coefficient in 2+1D""" - # Note: assumes pouch cell geometry + # Note: assumes pouch cell geometry return ( -2 * self.param.h / (self.param.delta ** 2) / self.param.l - 2 * (self.param.l_y + self.param.l_z) * self.param.h / self.param.delta diff --git a/tests/unit/test_models/test_submodels/test_interface/test_lead_acid.py b/tests/unit/test_models/test_submodels/test_interface/test_lead_acid.py index bad1837147..da98c47ac0 100644 --- a/tests/unit/test_models/test_submodels/test_interface/test_lead_acid.py +++ b/tests/unit/test_models/test_submodels/test_interface/test_lead_acid.py @@ -20,6 +20,7 @@ def test_public_functions(self): "Negative electrolyte potential": a_n, "Negative electrode open circuit potential": a_n, "Negative electrolyte concentration": a_n, + "Negative electrode temperature": a_n, } submodel = pybamm.interface.lead_acid.ButlerVolmer(param, "Negative") std_tests = tests.StandardSubModelTests(submodel, variables) @@ -32,6 +33,7 @@ def test_public_functions(self): "Positive electrolyte potential": a_p, "Positive electrode open circuit potential": a_p, "Positive electrolyte concentration": a_p, + "Positive electrode temperature": a_p, "Negative electrode interfacial current density": a_n, "Negative electrode exchange current density": a_n, } From 77778a6ff8cc1824f7dfbed5c8d7f0a8e53afb4c Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Wed, 30 Oct 2019 16:43:31 +0000 Subject: [PATCH 07/11] #678 fix add subtract simplify bug --- pybamm/expression_tree/binary_operators.py | 20 +++++++++++++++---- .../test_operations/test_simplify.py | 11 ++++++++++ 2 files changed, 27 insertions(+), 4 deletions(-) diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py index 70e9b1ea6d..0645d5b97b 100644 --- a/pybamm/expression_tree/binary_operators.py +++ b/pybamm/expression_tree/binary_operators.py @@ -277,9 +277,15 @@ def _binary_simplify(self, left, right): return left # Check matrices after checking scalars if is_matrix_zero(left): - return right + if isinstance(right, pybamm.Scalar): + return pybamm.Array(right.value * np.ones(left.shape_for_testing)) + else: + return right if is_matrix_zero(right): - return left + if isinstance(left, pybamm.Scalar): + return pybamm.Array(left.value * np.ones(right.shape_for_testing)) + else: + return left return pybamm.simplify_addition_subtraction(self.__class__, left, right) @@ -325,9 +331,15 @@ def _binary_simplify(self, left, right): return left # Check matrices after checking scalars if is_matrix_zero(left): - return -right + if isinstance(right, pybamm.Scalar): + return pybamm.Array(-right.value * np.ones(left.shape_for_testing)) + else: + return -right if is_matrix_zero(right): - return left + if isinstance(left, pybamm.Scalar): + return pybamm.Array(left.value * np.ones(right.shape_for_testing)) + else: + return left return pybamm.simplify_addition_subtraction(self.__class__, left, right) diff --git a/tests/unit/test_expression_tree/test_operations/test_simplify.py b/tests/unit/test_expression_tree/test_operations/test_simplify.py index 1d1f61c805..f84909706a 100644 --- a/tests/unit/test_expression_tree/test_operations/test_simplify.py +++ b/tests/unit/test_expression_tree/test_operations/test_simplify.py @@ -102,6 +102,17 @@ def myfunction(x, y): self.assertIsInstance((b - a).simplify(), pybamm.Scalar) self.assertEqual((b - a).simplify().evaluate(), 1) + # addition and subtraction with matrix zero + v = pybamm.Vector(np.zeros((10, 1))) + self.assertIsInstance((b + v).simplify(), pybamm.Array) + np.testing.assert_array_equal((b + v).simplify().evaluate(), np.ones((10, 1))) + self.assertIsInstance((v + b).simplify(), pybamm.Array) + np.testing.assert_array_equal((v + b).simplify().evaluate(), np.ones((10, 1))) + self.assertIsInstance((b - v).simplify(), pybamm.Array) + np.testing.assert_array_equal((b - v).simplify().evaluate(), np.ones((10, 1))) + self.assertIsInstance((v - b).simplify(), pybamm.Array) + np.testing.assert_array_equal((v - b).simplify().evaluate(), -np.ones((10, 1))) + # multiplication self.assertIsInstance((a * b).simplify(), pybamm.Scalar) self.assertEqual((a * b).simplify().evaluate(), 0) From 17b34e5b162b1c39a2b8efda0887b33c2c0f604b Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Wed, 30 Oct 2019 17:10:21 +0000 Subject: [PATCH 08/11] #687 resolving negative electrode heat discrepency --- .../compare_comsol_thermal.py | 158 ++++++++++++------ 1 file changed, 108 insertions(+), 50 deletions(-) diff --git a/results/comsol_comparison/compare_comsol_thermal.py b/results/comsol_comparison/compare_comsol_thermal.py index 9c47f8053b..66ba94f5dd 100644 --- a/results/comsol_comparison/compare_comsol_thermal.py +++ b/results/comsol_comparison/compare_comsol_thermal.py @@ -150,16 +150,10 @@ def myinterp(t): plt.legend() pybamm_voltage = pybamm.ProcessedVariable( - pybamm_model.variables["Terminal voltage [V]"], - solution.t, - solution.y, - mesh=mesh, + pybamm_model.variables["Terminal voltage [V]"], solution.t, solution.y, mesh=mesh )(plot_times / tau) comsol_voltage = pybamm.ProcessedVariable( - comsol_model.variables["Terminal voltage [V]"], - solution.t, - solution.y, - mesh=mesh, + comsol_model.variables["Terminal voltage [V]"], solution.t, solution.y, mesh=mesh )(plot_times / tau) plt.figure() plt.plot(plot_times, pybamm_voltage, "-", label="PyBaMM") @@ -175,39 +169,39 @@ def myinterp(t): x = mesh.combine_submeshes(*whole_cell)[0].nodes -def comparison_plot(var, plot_times=None): +def whole_cell_by_domain_comparison_plot(var, plot_times=None): """ - Plot pybamm heat source (defined over whole cell) against comsol heat source + Plot pybamm variable (defined over whole cell) against comsol variable (defined by component) """ if plot_times is None: plot_times = comsol_variables["time"] - # Process pybamm heat source - pybamm_q = pybamm.ProcessedVariable( + # Process pybamm variable + pybamm_var = pybamm.ProcessedVariable( pybamm_model.variables[var], solution.t, solution.y, mesh=mesh ) - # Process comsol heat source in negative electrode - comsol_q_n = pybamm.ProcessedVariable( + # Process comsol variable in negative electrode + comsol_var_n = pybamm.ProcessedVariable( comsol_model.variables["Negative electrode " + var[0].lower() + var[1:]], solution.t, solution.y, mesh=mesh, ) - # Process comsol heat source in separator (if defined here) + # Process comsol variable in separator (if defined here) try: - comsol_q_s = pybamm.ProcessedVariable( + comsol_var_s = pybamm.ProcessedVariable( comsol_model.variables["Separator " + var[0].lower() + var[1:]], solution.t, solution.y, mesh=mesh, ) except KeyError: - comsol_q_s = None + comsol_var_s = None print("Variable " + var + " not defined in separator") - # Process comsol heat source in positive electrode - comsol_q_p = pybamm.ProcessedVariable( + # Process comsol variable in positive electrode + comsol_var_p = pybamm.ProcessedVariable( comsol_model.variables["Positive electrode " + var[0].lower() + var[1:]], solution.t, solution.y, @@ -215,7 +209,7 @@ def comparison_plot(var, plot_times=None): ) # Make plot - if comsol_q_s: + if comsol_var_s: n_cols = 3 else: n_cols = 2 @@ -224,21 +218,21 @@ def comparison_plot(var, plot_times=None): for ind, t in enumerate(plot_times): color = cmap(float(ind) / len(plot_times)) - ax[0].plot(x_n * L_x, pybamm_q(x=x_n, t=t / tau), "-", color=color) - ax[0].plot(x_n * L_x, comsol_q_n(x=x_n, t=t / tau), "o", color=color) - if comsol_q_s: - ax[1].plot(x_s * L_x, pybamm_q(x=x_s, t=t / tau), "-", color=color) - ax[1].plot(x_s * L_x, comsol_q_s(x=x_s, t=t / tau), "o", color=color) + ax[0].plot(x_n * L_x, pybamm_var(x=x_n, t=t / tau), "-", color=color) + ax[0].plot(x_n * L_x, comsol_var_n(x=x_n, t=t / tau), "o", color=color) + if comsol_var_s: + ax[1].plot(x_s * L_x, pybamm_var(x=x_s, t=t / tau), "-", color=color) + ax[1].plot(x_s * L_x, comsol_var_s(x=x_s, t=t / tau), "o", color=color) ax[n_cols - 1].plot( x_p * L_x, - pybamm_q(x=x_p, t=t / tau), + pybamm_var(x=x_p, t=t / tau), "-", color=color, label="PyBaMM" if ind == 0 else "", ) ax[n_cols - 1].plot( x_p * L_x, - comsol_q_p(x=x_p, t=t / tau), + comsol_var_p(x=x_p, t=t / tau), "o", color=color, label="COMSOL" if ind == 0 else "", @@ -246,7 +240,7 @@ def comparison_plot(var, plot_times=None): ax[0].set_xlabel("x_n") ax[0].set_ylabel(var) - if comsol_q_s: + if comsol_var_s: ax[1].set_xlabel("x_s") ax[1].set_ylabel(var) ax[n_cols - 1].set_xlabel("x_p") @@ -255,28 +249,80 @@ def comparison_plot(var, plot_times=None): plt.tight_layout() -def temperature_plot(plot_times=None): +def electrode_comparison_plot(var, plot_times=None): """ - Plot pybamm heat source (defined over whole cell) against comsol heat source - (defined by component) + Plot pybamm variable against comsol variable (both defined separately in the + negative and positive electrode) """ if plot_times is None: plot_times = comsol_variables["time"] - # Process pybamm heat source - pybamm_T = pybamm.ProcessedVariable( - pybamm_model.variables["Cell temperature [K]"], - solution.t, - solution.y, - mesh=mesh, + # Process pybamm variable in negative electrode + pybamm_var_n = pybamm.ProcessedVariable( + pybamm_model.variables["Negative " + var], solution.t, solution.y, mesh=mesh ) - # Process comsol heat source in negative electrode - comsol_T = pybamm.ProcessedVariable( - comsol_model.variables["Cell temperature [K]"], - solution.t, - solution.y, - mesh=mesh, + # Process pybamm variable in positive electrode + pybamm_var_p = pybamm.ProcessedVariable( + pybamm_model.variables["Positive " + var], solution.t, solution.y, mesh=mesh + ) + + # Process comsol variable in negative electrode + comsol_var_n = pybamm.ProcessedVariable( + comsol_model.variables["Negative " + var], solution.t, solution.y, mesh=mesh + ) + + # Process comsol variable in positive electrode + comsol_var_p = pybamm.ProcessedVariable( + comsol_model.variables["Positive " + var], solution.t, solution.y, mesh=mesh + ) + + # Make plot + fig, ax = plt.subplots(1, 2, figsize=(15, 8)) + cmap = plt.get_cmap("inferno") + + for ind, t in enumerate(plot_times): + color = cmap(float(ind) / len(plot_times)) + ax[0].plot(x_n * L_x, pybamm_var_n(x=x_n, t=t / tau), "-", color=color) + ax[0].plot(x_n * L_x, comsol_var_n(x=x_n, t=t / tau), "o", color=color) + ax[1].plot( + x_p * L_x, + pybamm_var_p(x=x_p, t=t / tau), + "-", + color=color, + label="PyBaMM" if ind == 0 else "", + ) + ax[1].plot( + x_p * L_x, + comsol_var_p(x=x_p, t=t / tau), + "o", + color=color, + label="COMSOL" if ind == 0 else "", + ) + + ax[0].set_xlabel("x_n") + ax[0].set_ylabel(var) + ax[1].set_xlabel("x_p") + ax[1].set_ylabel(var) + plt.legend() + plt.tight_layout() + + +def whole_cell_comparison_plot(var, plot_times=None): + """ + Plot pybamm variable against comsol variable (both defined over whole cell) + """ + if plot_times is None: + plot_times = comsol_variables["time"] + + # Process pybamm variable + pybamm_var = pybamm.ProcessedVariable( + pybamm_model.variables[var], solution.t, solution.y, mesh=mesh + ) + + # Process comsol variable + comsol_var = pybamm.ProcessedVariable( + comsol_model.variables[var], solution.t, solution.y, mesh=mesh ) # Make plot @@ -287,29 +333,41 @@ def temperature_plot(plot_times=None): color = cmap(float(ind) / len(plot_times)) plt.plot( x * L_x, - pybamm_T(x=x, t=t / tau), + pybamm_var(x=x, t=t / tau), "-", color=color, label="PyBaMM" if ind == 0 else "", ) plt.plot( x * L_x, - comsol_T(x=x, t=t / tau), + comsol_var(x=x, t=t / tau), "o", color=color, label="COMSOL" if ind == 0 else "", ) plt.xlabel("x") - plt.ylabel("Temperature [K]") + plt.ylabel(var) plt.legend() plt.tight_layout() # Make plots plot_times = comsol_variables["time"][0::10] -comparison_plot("Irreversible electrochemical heating [W.m-3]", plot_times=plot_times) -comparison_plot("Reversible heating [W.m-3]", plot_times=plot_times) -comparison_plot("Total heating [W.m-3]", plot_times=plot_times) -# temperature_plot(plot_times) +# heat sources +whole_cell_by_domain_comparison_plot( + "Irreversible electrochemical heating [W.m-3]", plot_times=plot_times +) +whole_cell_by_domain_comparison_plot( + "Reversible heating [W.m-3]", plot_times=plot_times +) +whole_cell_by_domain_comparison_plot("Total heating [W.m-3]", plot_times=plot_times) +# potentials +electrode_comparison_plot("electrode potential [V]", plot_times=plot_times) +whole_cell_comparison_plot("Electrolyte potential [V]", plot_times=plot_times) +# concentrations +electrode_comparison_plot( + "particle surface concentration [mol.m-3]", plot_times=plot_times +) +whole_cell_comparison_plot("Electrolyte concentration [mol.m-3]", plot_times=plot_times) plt.show() From b377dadcdd7c74ba9f4cd04aa617077cc5c68c81 Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Thu, 31 Oct 2019 11:01:26 +0000 Subject: [PATCH 09/11] #678 add temp depend to backward tafel --- .../submodels/interface/kinetics/tafel.py | 4 +- .../compare_comsol_thermal.py | 373 ------------------ results/comsol_comparison/load_comsol_data.py | 142 ------- 3 files changed, 2 insertions(+), 517 deletions(-) delete mode 100644 results/comsol_comparison/compare_comsol_thermal.py delete mode 100644 results/comsol_comparison/load_comsol_data.py diff --git a/pybamm/models/submodels/interface/kinetics/tafel.py b/pybamm/models/submodels/interface/kinetics/tafel.py index a263657353..26ab18e807 100644 --- a/pybamm/models/submodels/interface/kinetics/tafel.py +++ b/pybamm/models/submodels/interface/kinetics/tafel.py @@ -82,5 +82,5 @@ class BackwardTafel(BaseModel): def __init__(self, param, domain): super().__init__(param, domain) - def _get_kinetics(self, j0, ne, eta_r): - return -j0 * pybamm.exp(-(ne / 2) * eta_r) + def _get_kinetics(self, j0, ne, eta_r, T): + return -j0 * pybamm.exp(-(ne / (2 * (1 + self.param.Theta * T))) * eta_r) diff --git a/results/comsol_comparison/compare_comsol_thermal.py b/results/comsol_comparison/compare_comsol_thermal.py deleted file mode 100644 index 66ba94f5dd..0000000000 --- a/results/comsol_comparison/compare_comsol_thermal.py +++ /dev/null @@ -1,373 +0,0 @@ -import pybamm -import numpy as np -import os -import pickle -import scipy.interpolate as interp -import matplotlib.pyplot as plt - -# change working directory to the root of pybamm -os.chdir(pybamm.root_dir()) - -"-----------------------------------------------------------------------------" -"Load comsol data" - -comsol_variables = pickle.load( - open("input/comsol_results/comsol_thermal_1C.pickle", "rb") -) - -"-----------------------------------------------------------------------------" -"Create and solve pybamm model" - -# load model and geometry -pybamm.set_logging_level("INFO") -options = {"thermal": "x-full"} -pybamm_model = pybamm.lithium_ion.DFN(options) -geometry = pybamm_model.default_geometry - -# load parameters and process model and geometry -param = pybamm_model.default_parameter_values -param.process_model(pybamm_model) -param.process_geometry(geometry) - -# create mesh -var = pybamm.standard_spatial_vars -var_pts = {var.x_n: 101, var.x_s: 51, var.x_p: 101, var.r_n: 31, var.r_p: 31} -mesh = pybamm.Mesh(geometry, pybamm_model.default_submesh_types, var_pts) - -# discretise model -disc = pybamm.Discretisation(mesh, pybamm_model.default_spatial_methods) -disc.process_model(pybamm_model) - -# discharge timescale -tau = param.process_symbol( - pybamm.standard_parameters_lithium_ion.tau_discharge -).evaluate(0) - -# solve model at comsol times -time = comsol_variables["time"] / tau -solver = pybamm.CasadiSolver() -solution = solver.solve(pybamm_model, time) - -"-----------------------------------------------------------------------------" -"Make Comsol 'model' for comparison" - -whole_cell = ["negative electrode", "separator", "positive electrode"] -comsol_t = comsol_variables["time"] -L_x = param.evaluate(pybamm.standard_parameters_lithium_ion.L_x) - - -def get_interp_fun(variable, domain): - """ - Create a :class:`pybamm.Function` object using the variable, to allow plotting with - :class:`'pybamm.QuickPlot'` (interpolate in space to match edges, and then create - function to interpolate in time) - """ - if domain == ["negative electrode"]: - comsol_x = comsol_variables["x_n"] - elif domain == ["separator"]: - comsol_x = comsol_variables["x_s"] - elif domain == ["positive electrode"]: - comsol_x = comsol_variables["x_p"] - elif domain == whole_cell: - comsol_x = comsol_variables["x"] - # Make sure to use dimensional space - pybamm_x = mesh.combine_submeshes(*domain)[0].nodes * L_x - variable = interp.interp1d(comsol_x, variable, axis=0)(pybamm_x) - - def myinterp(t): - return interp.interp1d(comsol_t, variable)(t)[:, np.newaxis] - - # Make sure to use dimensional time - fun = pybamm.Function(myinterp, pybamm.t * tau) - fun.domain = domain - return fun - - -comsol_c_n_surf = get_interp_fun(comsol_variables["c_n_surf"], ["negative electrode"]) -comsol_c_e = get_interp_fun(comsol_variables["c_e"], whole_cell) -comsol_c_p_surf = get_interp_fun(comsol_variables["c_p_surf"], ["positive electrode"]) -comsol_phi_n = get_interp_fun(comsol_variables["phi_n"], ["negative electrode"]) -comsol_phi_e = get_interp_fun(comsol_variables["phi_e"], whole_cell) -comsol_phi_p = get_interp_fun(comsol_variables["phi_p"], ["positive electrode"]) -comsol_voltage = interp.interp1d(comsol_t, comsol_variables["voltage"]) -comsol_temperature = get_interp_fun(comsol_variables["temperature"], whole_cell) -comsol_temperature_av = interp.interp1d( - comsol_t, comsol_variables["average temperature"] -) -comsol_q_irrev_n = get_interp_fun(comsol_variables["Q_irrev_n"], ["negative electrode"]) -comsol_q_irrev_p = get_interp_fun(comsol_variables["Q_irrev_p"], ["positive electrode"]) -comsol_q_rev_n = get_interp_fun(comsol_variables["Q_rev_n"], ["negative electrode"]) -comsol_q_rev_p = get_interp_fun(comsol_variables["Q_rev_p"], ["positive electrode"]) -comsol_q_total_n = get_interp_fun(comsol_variables["Q_total_n"], ["negative electrode"]) -comsol_q_total_s = get_interp_fun(comsol_variables["Q_total_s"], ["separator"]) -comsol_q_total_p = get_interp_fun(comsol_variables["Q_total_p"], ["positive electrode"]) - -# Create comsol model with dictionary of Matrix variables -comsol_model = pybamm.BaseModel() -comsol_model.variables = { - "Negative particle surface concentration [mol.m-3]": comsol_c_n_surf, - "Electrolyte concentration [mol.m-3]": comsol_c_e, - "Positive particle surface concentration [mol.m-3]": comsol_c_p_surf, - "Current [A]": pybamm_model.variables["Current [A]"], - "Negative electrode potential [V]": comsol_phi_n, - "Electrolyte potential [V]": comsol_phi_e, - "Positive electrode potential [V]": comsol_phi_p, - "Terminal voltage [V]": pybamm.Function(comsol_voltage, pybamm.t * tau), - "Cell temperature [K]": comsol_temperature, - "Volume-averaged cell temperature [K]": pybamm.Function( - comsol_temperature_av, pybamm.t * tau - ), - "Negative electrode irreversible electrochemical heating [W.m-3]": comsol_q_irrev_n, - "Positive electrode irreversible electrochemical heating [W.m-3]": comsol_q_irrev_p, - "Negative electrode reversible heating [W.m-3]": comsol_q_rev_n, - "Positive electrode reversible heating [W.m-3]": comsol_q_rev_p, - "Negative electrode total heating [W.m-3]": comsol_q_total_n, - "Separator total heating [W.m-3]": comsol_q_total_s, - "Positive electrode total heating [W.m-3]": comsol_q_total_p, -} - -"-----------------------------------------------------------------------------" -"Plot comparison" - -plot_times = comsol_variables["time"] -pybamm_T = pybamm.ProcessedVariable( - pybamm_model.variables["Volume-averaged cell temperature [K]"], - solution.t, - solution.y, - mesh=mesh, -)(plot_times / tau) -comsol_T = pybamm.ProcessedVariable( - comsol_model.variables["Volume-averaged cell temperature [K]"], - solution.t, - solution.y, - mesh=mesh, -)(plot_times / tau) -plt.figure() -plt.plot(plot_times, pybamm_T, "-", label="PyBaMM") -plt.plot(plot_times, comsol_T, "o", label="COMSOL") -plt.xlabel("t") -plt.ylabel("T") -plt.legend() - -pybamm_voltage = pybamm.ProcessedVariable( - pybamm_model.variables["Terminal voltage [V]"], solution.t, solution.y, mesh=mesh -)(plot_times / tau) -comsol_voltage = pybamm.ProcessedVariable( - comsol_model.variables["Terminal voltage [V]"], solution.t, solution.y, mesh=mesh -)(plot_times / tau) -plt.figure() -plt.plot(plot_times, pybamm_voltage, "-", label="PyBaMM") -plt.plot(plot_times, comsol_voltage, "o", label="COMSOL") -plt.xlabel("t") -plt.ylabel("Voltage [V]") -plt.legend() - -# Get mesh nodes -x_n = mesh.combine_submeshes(*["negative electrode"])[0].nodes -x_s = mesh.combine_submeshes(*["separator"])[0].nodes -x_p = mesh.combine_submeshes(*["positive electrode"])[0].nodes -x = mesh.combine_submeshes(*whole_cell)[0].nodes - - -def whole_cell_by_domain_comparison_plot(var, plot_times=None): - """ - Plot pybamm variable (defined over whole cell) against comsol variable - (defined by component) - """ - if plot_times is None: - plot_times = comsol_variables["time"] - - # Process pybamm variable - pybamm_var = pybamm.ProcessedVariable( - pybamm_model.variables[var], solution.t, solution.y, mesh=mesh - ) - - # Process comsol variable in negative electrode - comsol_var_n = pybamm.ProcessedVariable( - comsol_model.variables["Negative electrode " + var[0].lower() + var[1:]], - solution.t, - solution.y, - mesh=mesh, - ) - # Process comsol variable in separator (if defined here) - try: - comsol_var_s = pybamm.ProcessedVariable( - comsol_model.variables["Separator " + var[0].lower() + var[1:]], - solution.t, - solution.y, - mesh=mesh, - ) - except KeyError: - comsol_var_s = None - print("Variable " + var + " not defined in separator") - # Process comsol variable in positive electrode - comsol_var_p = pybamm.ProcessedVariable( - comsol_model.variables["Positive electrode " + var[0].lower() + var[1:]], - solution.t, - solution.y, - mesh=mesh, - ) - - # Make plot - if comsol_var_s: - n_cols = 3 - else: - n_cols = 2 - fig, ax = plt.subplots(1, n_cols, figsize=(15, 8)) - cmap = plt.get_cmap("inferno") - - for ind, t in enumerate(plot_times): - color = cmap(float(ind) / len(plot_times)) - ax[0].plot(x_n * L_x, pybamm_var(x=x_n, t=t / tau), "-", color=color) - ax[0].plot(x_n * L_x, comsol_var_n(x=x_n, t=t / tau), "o", color=color) - if comsol_var_s: - ax[1].plot(x_s * L_x, pybamm_var(x=x_s, t=t / tau), "-", color=color) - ax[1].plot(x_s * L_x, comsol_var_s(x=x_s, t=t / tau), "o", color=color) - ax[n_cols - 1].plot( - x_p * L_x, - pybamm_var(x=x_p, t=t / tau), - "-", - color=color, - label="PyBaMM" if ind == 0 else "", - ) - ax[n_cols - 1].plot( - x_p * L_x, - comsol_var_p(x=x_p, t=t / tau), - "o", - color=color, - label="COMSOL" if ind == 0 else "", - ) - - ax[0].set_xlabel("x_n") - ax[0].set_ylabel(var) - if comsol_var_s: - ax[1].set_xlabel("x_s") - ax[1].set_ylabel(var) - ax[n_cols - 1].set_xlabel("x_p") - ax[n_cols - 1].set_ylabel(var) - plt.legend() - plt.tight_layout() - - -def electrode_comparison_plot(var, plot_times=None): - """ - Plot pybamm variable against comsol variable (both defined separately in the - negative and positive electrode) - """ - if plot_times is None: - plot_times = comsol_variables["time"] - - # Process pybamm variable in negative electrode - pybamm_var_n = pybamm.ProcessedVariable( - pybamm_model.variables["Negative " + var], solution.t, solution.y, mesh=mesh - ) - - # Process pybamm variable in positive electrode - pybamm_var_p = pybamm.ProcessedVariable( - pybamm_model.variables["Positive " + var], solution.t, solution.y, mesh=mesh - ) - - # Process comsol variable in negative electrode - comsol_var_n = pybamm.ProcessedVariable( - comsol_model.variables["Negative " + var], solution.t, solution.y, mesh=mesh - ) - - # Process comsol variable in positive electrode - comsol_var_p = pybamm.ProcessedVariable( - comsol_model.variables["Positive " + var], solution.t, solution.y, mesh=mesh - ) - - # Make plot - fig, ax = plt.subplots(1, 2, figsize=(15, 8)) - cmap = plt.get_cmap("inferno") - - for ind, t in enumerate(plot_times): - color = cmap(float(ind) / len(plot_times)) - ax[0].plot(x_n * L_x, pybamm_var_n(x=x_n, t=t / tau), "-", color=color) - ax[0].plot(x_n * L_x, comsol_var_n(x=x_n, t=t / tau), "o", color=color) - ax[1].plot( - x_p * L_x, - pybamm_var_p(x=x_p, t=t / tau), - "-", - color=color, - label="PyBaMM" if ind == 0 else "", - ) - ax[1].plot( - x_p * L_x, - comsol_var_p(x=x_p, t=t / tau), - "o", - color=color, - label="COMSOL" if ind == 0 else "", - ) - - ax[0].set_xlabel("x_n") - ax[0].set_ylabel(var) - ax[1].set_xlabel("x_p") - ax[1].set_ylabel(var) - plt.legend() - plt.tight_layout() - - -def whole_cell_comparison_plot(var, plot_times=None): - """ - Plot pybamm variable against comsol variable (both defined over whole cell) - """ - if plot_times is None: - plot_times = comsol_variables["time"] - - # Process pybamm variable - pybamm_var = pybamm.ProcessedVariable( - pybamm_model.variables[var], solution.t, solution.y, mesh=mesh - ) - - # Process comsol variable - comsol_var = pybamm.ProcessedVariable( - comsol_model.variables[var], solution.t, solution.y, mesh=mesh - ) - - # Make plot - plt.figure(figsize=(15, 8)) - cmap = plt.get_cmap("inferno") - - for ind, t in enumerate(plot_times): - color = cmap(float(ind) / len(plot_times)) - plt.plot( - x * L_x, - pybamm_var(x=x, t=t / tau), - "-", - color=color, - label="PyBaMM" if ind == 0 else "", - ) - plt.plot( - x * L_x, - comsol_var(x=x, t=t / tau), - "o", - color=color, - label="COMSOL" if ind == 0 else "", - ) - - plt.xlabel("x") - plt.ylabel(var) - plt.legend() - plt.tight_layout() - - -# Make plots -plot_times = comsol_variables["time"][0::10] -# heat sources -whole_cell_by_domain_comparison_plot( - "Irreversible electrochemical heating [W.m-3]", plot_times=plot_times -) -whole_cell_by_domain_comparison_plot( - "Reversible heating [W.m-3]", plot_times=plot_times -) -whole_cell_by_domain_comparison_plot("Total heating [W.m-3]", plot_times=plot_times) -# potentials -electrode_comparison_plot("electrode potential [V]", plot_times=plot_times) -whole_cell_comparison_plot("Electrolyte potential [V]", plot_times=plot_times) -# concentrations -electrode_comparison_plot( - "particle surface concentration [mol.m-3]", plot_times=plot_times -) -whole_cell_comparison_plot("Electrolyte concentration [mol.m-3]", plot_times=plot_times) -plt.show() diff --git a/results/comsol_comparison/load_comsol_data.py b/results/comsol_comparison/load_comsol_data.py deleted file mode 100644 index a67ab9b8f3..0000000000 --- a/results/comsol_comparison/load_comsol_data.py +++ /dev/null @@ -1,142 +0,0 @@ -import pybamm -import os -import pandas as pd -import pickle -import numpy as np - -# change working directory the root of pybamm -os.chdir(pybamm.root_dir()) - -# set filepath for data and name of file to pickle to -path = "input/comsol_results_csv/thermal/1C/" -savefile = "input/comsol_results/comsol_thermal_1C.pickle" - -# time-voltage -comsol = pd.read_csv(path + "voltage.csv", sep=",", header=None) -comsol_time = comsol[0].values -comsol_time_npts = len(comsol_time) -comsol_voltage = comsol[1].values - -# negative electrode potential -comsol = pd.read_csv(path + "phi_s_n.csv", sep=",", header=None) -comsol_x_n_npts = int(len(comsol[0].values) / comsol_time_npts) -comsol_x_n = comsol[0].values[0:comsol_x_n_npts] -comsol_phi_n_vals = np.reshape( - comsol[1].values, (comsol_x_n_npts, comsol_time_npts), order="F" -) - -# negative particle surface concentration -comsol = pd.read_csv(path + "c_n_surf.csv", sep=",", header=None) -comsol_c_n_surf_vals = np.reshape( - comsol[1].values, (comsol_x_n_npts, comsol_time_npts), order="F" -) - -# positive electrode potential -comsol = pd.read_csv(path + "phi_s_p.csv", sep=",", header=None) -comsol_x_p_npts = int(len(comsol[0].values) / comsol_time_npts) -comsol_x_p = comsol[0].values[0:comsol_x_p_npts] -comsol_phi_p_vals = np.reshape( - comsol[1].values, (comsol_x_p_npts, comsol_time_npts), order="F" -) - -# positive particle surface concentration -comsol = pd.read_csv(path + "c_p_surf.csv", sep=",", header=None) -comsol_c_p_surf_vals = np.reshape( - comsol[1].values, (comsol_x_p_npts, comsol_time_npts), order="F" -) - -# electrolyte concentration -comsol = pd.read_csv(path + "c_e.csv", sep=",", header=None) -comsol_x_npts = int(len(comsol[0].values) / comsol_time_npts) -comsol_x = comsol[0].values[0:comsol_x_npts] -comsol_c_e_vals = np.reshape( - comsol[1].values, (comsol_x_npts, comsol_time_npts), order="F" -) - -# electrolyte potential -comsol = pd.read_csv(path + "phi_e.csv", sep=",", header=None) -comsol_phi_e_vals = np.reshape( - comsol[1].values, (comsol_x_npts, comsol_time_npts), order="F" -) - -# temperature -comsol = pd.read_csv(path + "temperature.csv", sep=",", header=None) -comsol_temp_vals = np.reshape( - comsol[1].values, (comsol_x_npts, comsol_time_npts), order="F" -) - -# average temperature -comsol_temp_av = np.mean(comsol_temp_vals, axis=0) - -# irreversible heat source in negative electrode -comsol = pd.read_csv(path + "q_irrev_n.csv", sep=",", header=None) -comsol_q_irrev_n_vals = np.reshape( - comsol[1].values, (comsol_x_n_npts, comsol_time_npts), order="F" -) - -# irreversible heat source in positive electrode -comsol = pd.read_csv(path + "q_irrev_p.csv", sep=",", header=None) -comsol_q_irrev_p_vals = np.reshape( - comsol[1].values, (comsol_x_p_npts, comsol_time_npts), order="F" -) - -# reversible heat source in negative electrode -comsol = pd.read_csv(path + "q_rev_n.csv", sep=",", header=None) -comsol_q_rev_n_vals = np.reshape( - comsol[1].values, (comsol_x_n_npts, comsol_time_npts), order="F" -) - -# reversible heat source in positive electrode -comsol = pd.read_csv(path + "q_rev_p.csv", sep=",", header=None) -comsol_q_rev_p_vals = np.reshape( - comsol[1].values, (comsol_x_p_npts, comsol_time_npts), order="F" -) - -# total heat source in negative electrode -comsol = pd.read_csv(path + "q_total_n.csv", sep=",", header=None) -comsol_q_total_n_vals = np.reshape( - comsol[1].values, (comsol_x_n_npts, comsol_time_npts), order="F" -) - -# total heat source in separator -comsol = pd.read_csv(path + "q_total_s.csv", sep=",", header=None) -comsol_x_s_npts = int(len(comsol[0].values) / comsol_time_npts) -comsol_x_s = comsol[0].values[0:comsol_x_s_npts] -comsol_q_total_s_vals = np.reshape( - comsol[1].values, (comsol_x_s_npts, comsol_time_npts), order="F" -) - -# total heat source in positive electrode -comsol = pd.read_csv(path + "q_total_p.csv", sep=",", header=None) -comsol_q_total_p_vals = np.reshape( - comsol[1].values, (comsol_x_p_npts, comsol_time_npts), order="F" -) - - -# add comsol variables to dict and pickle -comsol_variables = { - "time": comsol_time, - "x_n": comsol_x_n, - "x_s": comsol_x_s, - "x_p": comsol_x_p, - "x": comsol_x, - "voltage": comsol_voltage, - "phi_n": comsol_phi_n_vals, - "phi_p": comsol_phi_p_vals, - "phi_e": comsol_phi_e_vals, - "c_n_surf": comsol_c_n_surf_vals, - "c_p_surf": comsol_c_p_surf_vals, - "c_e": comsol_c_e_vals, - "temperature": comsol_temp_vals, - "average temperature": comsol_temp_av, - "Q_irrev_n": comsol_q_irrev_n_vals, - "Q_irrev_p": comsol_q_irrev_p_vals, - "Q_rev_n": comsol_q_rev_n_vals, - "Q_rev_p": comsol_q_rev_p_vals, - "Q_total_n": comsol_q_total_n_vals, - "Q_total_s": comsol_q_total_s_vals, - "Q_total_p": comsol_q_total_p_vals, -} - -with open(savefile, "wb") as f: - pickle.dump(comsol_variables, f, pickle.HIGHEST_PROTOCOL) From b89df67144192fc6f28e727d87a4cc1e4b4c9a5b Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Thu, 31 Oct 2019 11:06:39 +0000 Subject: [PATCH 10/11] #678 update changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 80b9c20626..52d195e7a4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,6 +19,7 @@ ## Bug fixes +- Adds missing temperature dependence in electrolyte and interface submodels ([#698](https://github.com/pybamm-team/PyBaMM/pull/698)) - Fix differentiation of functions that have more than one argument ([#687](https://github.com/pybamm-team/PyBaMM/pull/687)) - Add warning if `ProcessedVariable` is called outisde its interpolation range ([#681](https://github.com/pybamm-team/PyBaMM/pull/681)) - Improve the way `ProcessedVariable` objects are created in higher dimensions ([#581](https://github.com/pybamm-team/PyBaMM/pull/581)) From 379c957fd2b6b59a5ee25a5e70697ea5e38b41e6 Mon Sep 17 00:00:00 2001 From: Robert Timms Date: Thu, 31 Oct 2019 12:29:01 +0000 Subject: [PATCH 11/11] #678 fix set temperature script --- results/2plus1D/set_temperature_spm_1plus1D.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/results/2plus1D/set_temperature_spm_1plus1D.py b/results/2plus1D/set_temperature_spm_1plus1D.py index fca61facbd..f6101a8824 100644 --- a/results/2plus1D/set_temperature_spm_1plus1D.py +++ b/results/2plus1D/set_temperature_spm_1plus1D.py @@ -66,7 +66,7 @@ def non_dim_temperature(temperature): "Negative current collector potential [V]", "Positive current collector potential [V]", "Current [A]", - "X-averaged total heating [A.V.m-3]", + "X-averaged total heating [W.m-3]", "X-averaged positive particle surface concentration [mol.m-3]", "X-averaged cell temperature [K]", ] @@ -145,16 +145,16 @@ def non_dim_temperature(temperature): for bat_id in range(nbat): plt.plot( solution1.t * t_hour, - processed_vars_step1["X-averaged total heating [A.V.m-3]"](solution1.t, z=z)[ + processed_vars_step1["X-averaged total heating [W.m-3]"](solution1.t, z=z)[ bat_id, : ], solution2.t * t_hour, - processed_vars_step2["X-averaged total heating [A.V.m-3]"](solution2.t, z=z)[ + processed_vars_step2["X-averaged total heating [W.m-3]"](solution2.t, z=z)[ bat_id, : ], ) plt.xlabel("t [hrs]") -plt.ylabel("X-averaged total heating [A.V.m-3]") +plt.ylabel("X-averaged total heating [W.m-3]") plt.yscale("log") # local concentration