Updating temperature in solver.step() - Setting up external thermal model

[ChaddKiggins]

Hi all,

I am trying to set up a thermal model where I can calculate the temperature gradient along the radial dimension of a cylindrical cell. I am currently making quite a few assumptions, but essentially I want to run a P2D at each wind along the radial dimension. The P2D's would provide the heat source and a simple thermal model would update the temperature at each wind position. I'm getting stuck on the next part where I take this updated temperature and plug it back into the next time step.

The code is below. The while loop is where the P2D's and thermal model is being solved. What I'm trying to do is get the T array (temperature after the thermal model) to feed back into the next solver.step(). I'm not to sure how to update the solution object to incorporate this temperature feedback.

Does anyone have any thoughts on how I can approach this?

import pybamm
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd


#Class to instantiate any given number of P2D models 
class P2dnode:
    def __init__(self, parameters):
        self.model = pybamm.lithium_ion.DFN({"thermal":"lumped"})
        self.parameters = pybamm.ParameterValues(f"{parameters}")
        self.step_solution = None
        
    
    def setup_model(self):
        self.parameters.process_model(self.model)

        self.geometry = self.model.default_geometry
        self.parameters.process_geometry(self.geometry)

        self.var = pybamm.standard_spatial_vars
        self.var_pts = {self.var.x_n: 20, self.var.x_s: 15, self.var.x_p: 20, self.var.r_n: 30, self.var.r_p: 30}
        self.mesh = pybamm.Mesh(self.geometry, self.model.default_submesh_types, self.var_pts)

        self.disc = pybamm.Discretisation(self.mesh, self.model.default_spatial_methods)
        self.disc.process_model(self.model)

    
#Instantiate a given number of P2D models
p2ds = [] #List to hold all the P2D models
for _ in range(10):
    model = P2dnode(parameters = "Chen2020")
    model.parameters.update({"Current function [A]": "[input]", "Volume-averaged cell temperature [K]": "[input]"}, check_already_exists=False)
    model.setup_model()
    p2ds.append(model)


#Setup mesh for the thermal model
L = 10e-3 #Length
N = 10 #Number of nodes (equal to the number of P2D's)
dr = L/N #Grid Size
r = np.linspace(dr/2, L-dr/2, N) #radius grid
dtt = 0.001 #Time step
t_final = 10 #Final time
Tin = 298.15 #INitial temp
rho = 2682 #Density
cp = 800 #Heat capacity
k = 0.5 #Thermal conductivity
h = 20 #heat transfer coefficient
T_ambient = 298.15 #ambient temperature        
alpha = k / (rho * cp) #Thermal diffusivity

T = np.ones(N) * Tin #Intial temperature
t = np.arange(0, t_final, dtt) #Time grid


#Setup Time intervals for stepping the solution
ti = 0
dt = 10
tf = 100

#While loop to timestep the solutions
while ti < tf:
    #Step 1: Solve P2D models
    Q = []
    for i, p2d in enumerate(p2ds):
        solver = pybamm.CasadiSolver()
        inputs = {"Current function [A]": 10, "Volume-averaged cell temperature [K]": T[i]}
        p2d.step_solution = solver.step(model=p2d.model, old_solution=p2d.step_solution, dt=dt, inputs=inputs)
        print(p2d.step_solution["Volume-averaged cell temperature [K]"].data[-1])
        q = p2d.step_solution["Total heating [W.m-3]"].data
        q = np.average(q)
        Q.append(q)
    print(Q)
    #Step 2: Solve 1D thermal model
    for j in range(1, len(t)):
        T_new = T.copy()
        T_new[0] = T[1] #Set up no heat transfer boundary condition
        for i in range(1, N-1):
            T_new[i] = T[i] + dtt * alpha * (((T[i+1] - 2*T[i] + T[i-1]) / dr**2 +  (T[i+1] - T[i-1]) / (2*r[i]*dr)) + Q[i]/k) #Discretize 1D heat equation in cylindrical coordiantes (radial)
        T_new[-1] = (k * T[-2] + h * dr * T_ambient) / (k + h *dr) #External boundary condition (convection)
        T = T_new #Update temperature for the next time step (feed back into the next round of P2D's)
    print(T)
    ti += dt
    

r = r * 1000 #Conver to millimeter
T = T - 273.15 #Convert back to Celcius

plt.plot(r, T)
plt.ylabel("Temperature (C)")
plt.xlabel("Radial Distance from Center (mm)")
plt.title("1D Radial Temp Generation")

plt.show()

[valentinsulzer]

"Volume-averaged cell temperature [K]" isn't a parameter, since you used check_already_exists=False it won't error but it will have no effect. You could try updating the ambient temperature instead (each cell responds to the average ambient temperature of its neighbors?), or do something with model.set_initial_conditions_from to update the initial temperature at each step, or use an isothermal model with {"calculate heat source for isothermal models": "true"} and then update the ambient temperature at each step (for an isothermal model the temperature is equal to the ambient temperature) and use the heat fluxes to calculate the heating in your 1D thermal model

@TomTranter is the real expert here so I'll let him weigh in

[ChaddKiggins]

Thanks @tinosulzer for the suggestions. Yeah that's what I figured with the "Volume-averaged cell temperature [K]", wasn't having any influence.

I like the idea of the model.set_initial_conditions_from. I'll see if I can reset the temperatures at the beginning of each solver.step() while retaining the rest of the Solution object.

Using the ambient temperature is very interesting. I'm not updating any temperatures directly from the P2D model's, just getting the heat source term. So updating the temps externally and then updating the ambient temps could work quite well if that can be read into the next step as the starting point. Realistically, all I care about is the temp that the next P2D starts on so I'll for sure give this a try.

[TomTranter]

Hey @ChaddKiggins this approach is the one I use in the jellyroll papers 10.1149/1945-7111/aba44b 10.1149/1945-7111/abd44f and I have the code for it here https://github.com/TomTranter/pybamm_pnm I coupled the PyBaMM models with a network solver called OpenPNM for solving the current and thermal around the rolls and updated the temperatures then stepped the electrochemical sims in parallel and updated internal resistances etc and heat sources. Solving back and forth it runs a single discharge with 1000 models in about 10 minutes. I then took the same logic and put it into a pack simulation which became liionpack. https://liionpack.readthedocs.io/en/latest/ You could potentially use liionpack without too much alteration to your code and might save you some work getting parallel processing etc up and running if you decide you need it. I think the relevant bit you could copy and paste into your simulation is here https://github.com/pybamm-team/liionpack/blob/cb38ed5113452d82ac4a63ec3ef34e62f61ad8ed/liionpack/simulations.py#L90
If you get stuck and need some help then message me on our slack channel and I'll try to help. The full Jellyroll code needs a bit of tidying up - however might serve as a useful comparison for you and I could help you get set up with it

[ChaddKiggins]

Thanks @TomTranter, was able to quickly modify my code to get the temperature from the previous time step incorporated into the next round of P2D's. Below is the code and some plots. With the temperature being different, the discharge at each node point is indeed different which is what I was aiming for (hard to see due to the small temp delta).

Thanks also for sending over all the other links. The Jellyroll code looks like a really great way to handle the local P2D's and apply whatever global model we may want, enabling an accurate 3D depiction of the cell. Would love to try it out. I'll take a deeper look into that and liionpack as well. After familiarizing myself with those, I'll reach out if I have questions (I'm sure I will) so I appreciate the help!

import pybamm
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib.colors import Normalize
from matplotlib.patches import Circle
from matplotlib import cm
from SDI_40T3_params import get_parameter_values_40T3_ChargeOpt



#Class to instantiate any given number of P2D models 
class P2dnode:
    def __init__(self, parameters):
        self.model = pybamm.lithium_ion.DFN({"thermal":"isothermal", "calculate heat source for isothermal models": "true"})
        self.parameters = pybamm.ParameterValues(parameters)
        self.parameters["Ambient temperature [K]"] = pybamm.InputParameter("Input temperature [K]")
        self.parameters["Initial temperature [K]"] = pybamm.InputParameter("Input temperature [K]")
        self.step_solution = None
        
    
    def setup_model(self):
        self.parameters.process_model(self.model)

        self.geometry = self.model.default_geometry
        self.parameters.process_geometry(self.geometry)

        self.var = pybamm.standard_spatial_vars
        self.var_pts = {self.var.x_n: 20, self.var.x_s: 15, self.var.x_p: 20, self.var.r_n: 30, self.var.r_p: 30}
        self.mesh = pybamm.Mesh(self.geometry, self.model.default_submesh_types, self.var_pts)

        self.disc = pybamm.Discretisation(self.mesh, self.model.default_spatial_methods)
        self.disc.process_model(self.model)

    
#Instantiate a given number of P2D models
p2ds = [] #List to hold all the P2D models
for _ in range(20):
    model = P2dnode(parameters = get_parameter_values_40T3_ChargeOpt())
    model.parameters.update({"Current function [A]": "[input]", "Ambient temperature [K]": "[input]", "Initial temperature [K]":"[input]"})
    model.parameters.set_initial_stoichiometries(0.99)
    model.setup_model()
    p2ds.append(model)


#Setup mesh for the thermal model
L = 10e-3 #Length
N = 20 #Number of nodes (equal to the number of P2D's)
dr = L/N #Grid Size
r = np.linspace(dr/2, L-dr/2, N) #radius grid
dtt = 0.01 #Time step
t_final = 10 #Final time
Tin = 298.15 #INitial temp
rho = 2682 #Density
cp = 800 #Heat capacity
k = 0.5 #Thermal conductivity
h = 10 #heat transfer coefficient
T_ambient = 298.15 #ambient temperature        
alpha = k / (rho * cp) #Thermal diffusivity

T = np.ones(N) * Tin #Intial temperature
t = np.arange(0, t_final, dtt) #Time grid


#Setup Time intervals for stepping the solution
ti = 0
dt = 10
tf = 720

#While loop to timestep the solutions
T_total = []
while ti < tf:
    #Step 1: Solve P2D models
    Q = []
    for i, p2d in enumerate(p2ds):
        solver = pybamm.CasadiSolver(mode="safe", atol=1e-06, rtol=1e-06, dt_max=10)
        inputs = {"Current function [A]": 20, "Ambient temperature [K]": T[i], "Initial temperature [K]": T[i]}
        p2d.step_solution = solver.step(model=p2d.model, old_solution=p2d.step_solution, dt=dt, inputs=inputs, save=True)
        print(p2d.step_solution["Ambient temperature [K]"].data[-1])
        q = p2d.step_solution["Total heating [W.m-3]"].data
        q = np.average(q)
        Q.append(q)
    print(Q)
    #Step 2: Solve 1D thermal model
    for j in range(1, len(t)):
        T_new = T.copy()
        T_new[0] = T[1] #Set up no heat transfer boundary condition (placeholder for now)
        for i in range(1, N-1):
            T_new[i] = T[i] + dtt * alpha * (((T[i+1] - 2*T[i] + T[i-1]) / dr**2 +  (T[i+1] - T[i-1]) / (2*r[i]*dr)) + Q[i]/k) #Discretize 1D heat equation in cylindrical coordiantes (radial)
        T_new[-1] = (k * T[-2] + h * dr * T_ambient) / (k + h *dr) #External boundary condition (convection)
        T = T_new #Update temperature for the next time step (feed back into the next round of P2D's)
    print(T)
    ti += dt

#df = pd.DataFrame(T_total)
#print(df)

r = r * 1000 #Conver to millimeter
T = T - 273.15 #Convert back to Celcius

plt.plot(r, T)
plt.ylabel("Temperature (C)")
plt.xlabel("Radial Distance from Center (mm)")
plt.title("1D Radial Temp Generation")


min = np.min(T)
max = np.max(T)

cmap = cm.get_cmap("rainbow")
norm = Normalize(vmin = min, vmax = max)

#Draw circle representing temp gradient in cylindrical cell

fig, ax = plt.subplots()


for i, rad in reversed(list(enumerate(r))):
    r_norm = r[i]/r[-1]
    color = cmap(norm(T[i]))
    circle = Circle((0.5,0.5), r_norm/2, color=color)
    ax.add_patch(circle)
    ax.axis("equal")
    ax.axis("off")

#circle = Circle((0.5, 0.5), r[-1]/r[-1], color = cmap(norm(T[15]))) #single circle for illustrative purposes (averaged case)
#ax.add_patch(circle)
#ax.axis("equal")
#ax.axis("off")

#Add core of cylindrical cell
circle_core = Circle((0.5,0.5), (r[1]/r[-1]/3), color="white")
ax.add_patch(circle_core)

#Add Colorbar for temp distribution
cbar = plt.colorbar(plt.cm.ScalarMappable(norm = norm, cmap=cmap), ax=ax)
cbar.ax.set_title("Temperature (C)")

plt.figure()
#check voltage curves with temperature feedback
vol = []
time = []
for p2d in p2ds:
    V = p2d.step_solution["Terminal voltage [V]"].data
    vol.append(V)
    t = p2d.step_solution["Time [s]"].data
    time.append(t)

plt.plot(time[0], vol[0], label="core")
plt.plot(time[-1], vol[-1], label = "case")
plt.ylabel("Voltage (V)")
plt.xlabel("Time (s)")
plt.title("Discharge Between Core and Case")

plt.legend()

plt.show()