How to set interface boundary conditions?

[xyyrly]

Hi,
I am creating a thermal model considering the temperature difference between active particle temperature and electrolyte. I hope to set a BC at particle surface like:

I tried to write like this (ignore heat sorce first):

But it returned an error, said"ModelError: model is underdetermined (too many variables)"
Are there some ways I can try to realize it? Thanks a loooot!

[brosaplanella]

As the error says, you have too many variables and too few equations so the system is underdetermined. It is a bit hard to miss which equation is missing, it would be great if you can share a minimum working example.

[xyyrly]

Thank you, brosaplanella!!!! I have fixed that problem, But there is another one.. as you see, the BC at surface is related to the heat flux difference, but the two heat flux iare in r direction and x direction respectively. I set the bc in pybamm by:

# boundary conditions 
T1 = pybamm.Parameter("T1 [K]")
T2 = pybamm.Parameter("T2 [K]")
grad_Tb_cathode_particle_right=lumbda_n_eff/lumbda_s*pybamm.grad(T_e_p)
model.boundary_conditions = {
    T_e: {"left": (T1, "Dirichlet"), "right": (T2, "Dirichlet")},
    Ts: {"left": (pybamm.Scalar(0), "Neumann"), "right": (grad_Tb_cathode_particle_right, "Neumann")}
}

# initial conditions
model.initial_conditions = {
    T_e: T1,
    Ts: T1
}

#Error
pybamm gave an error, I think it is beacuse the two heat flux vectors are not in the same direction. Are there some way I can change pybamm.grad to its norm?

Overwhelmed with gratitude~~

[brosaplanella]

Can you please share a minimum working example (as in a full script that has the minimum code to reproduce the error) so we can run it at our end?

[xyyrly]

Yes, Thanks a lot!!

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

model = pybamm.BaseModel()
T_e_s = pybamm.Variable("Separator electrolyte temperature [K]", domain="separator")
T_e_p = pybamm.Variable("Positive electrolyte temperature [K]", domain="positive electrode")
T_e = pybamm.concatenation(T_e_s, T_e_p)
Ts = pybamm.Variable(
    "Positive particle temperature [K]",
    domain="positive particle",
    auxiliary_domains={
        "secondary": "positive electrode",
    },
)
from scipy import constants
L_s = pybamm.Parameter("Separator thickness [m]")
L_p = pybamm.Parameter("Positive electrode thickness [m]")
eps_n = pybamm.PrimaryBroadcast(
    pybamm.Parameter("Positive electrode porosity [1]"), "positive electrode"
)
eps_sep = pybamm.PrimaryBroadcast(
    pybamm.Parameter("Separator porosity [1]"), "separator"
)
eps = pybamm.concatenation(eps_sep, eps_n)

rho_n = pybamm.PrimaryBroadcast(
    pybamm.Parameter("Positive electrode density [kg.m-3]"), "positive electrode"
)
rho_sep = pybamm.PrimaryBroadcast(
    pybamm.Parameter("Separator density [kg.m-3]"), "separator"
)
rho = pybamm.concatenation(rho_sep, rho_n)

lumbda_n = pybamm.PrimaryBroadcast(
    pybamm.Parameter("Positive electrode heat conductivity [W.m-2.K-1]"), "positive electrode"
)
lumbda_sep = pybamm.PrimaryBroadcast(
    pybamm.Parameter("Separator heat conductivity [W.m-2.K-1]"), "separator"
)
lumbda = pybamm.concatenation(lumbda_sep, lumbda_n)

cp_n = pybamm.PrimaryBroadcast(
    pybamm.Parameter("Positive electrode heat capacity [J.kg-1.K-1]"), "positive electrode"
)
cp_sep = pybamm.PrimaryBroadcast(
    pybamm.Parameter("Separator heat capacity [J.kg-1.K-1]"), "separator"
)
cp = pybamm.concatenation(cp_sep, cp_n)

R_p = pybamm.Parameter("Positive particle radius [m]")
rho_s = pybamm.Parameter("Positive particle density [kg.m-3]")
cp_s = pybamm.Parameter("Positive particle heat capacity [J.kg-1.K-1]")
lumbda_s = pybamm.Parameter("Positive particle heat conductivity [W.m-2.K-1]")
# eqn
alpha_n=1/eps_n/rho_n/cp_n
alpha_sep=1/rho_sep/cp_sep
alpha = pybamm.concatenation(alpha_sep, alpha_n)
lumbda_n_eff=lumbda_n*(eps_n**1.5)
lumbda_sep_eff=lumbda_sep
lumbda = pybamm.concatenation(lumbda_sep, lumbda_n)

N=lumbda*pybamm.grad(T_e)
dTdt=alpha*pybamm.div(N)
N_particle=lumbda_s*pybamm.grad(Ts)
dTsdt=pybamm.div(N_particle)/rho_s/cp_s
model.rhs={T_e:dTdt,
           Ts:dTsdt,
          }
T1 = pybamm.Parameter("T1 [K]")
T2 = pybamm.Parameter("T2 [K]")
grad_Tb_cathode_particle_right=lumbda_n_eff/lumbda_s*pybamm.grad(T_e_p)
# boundary conditions 
model.boundary_conditions = {
    T_e: {"left": (T1, "Dirichlet"), "right": (T2, "Dirichlet")},
    Ts: {"left": (pybamm.Scalar(0), "Neumann"), "right": (grad_Tb_cathode_particle_right, "Neumann")}
}

# initial conditions
model.initial_conditions = {
    T_e: T1,
    Ts: T1
}
model.variables={
    "electrolyte temperature [K]": T_e,
    "Positive particle temperature [K]": Ts,
}
param = pybamm.ParameterValues(
    {
        "T1 [K]":300,
        "T2 [K]":400,
        "T3 [K]":300,
        "Separator thickness [m]": 25e-6,
        "Positive electrode thickness [m]": 100e-6,
        "Positive electrode porosity [1]": 0.417,
        "Positive electrode density [kg.m-3]": 1000,
        "Positive electrode heat capacity [J.kg-1.K-1]": 700,
        "Positive electrode heat conductivity [W.m-2.K-1]": 0.16,
        "Separator porosity [1]":0.417,
        "Separator density [kg.m-3]":1100,
        "Separator heat capacity [J.kg-1.K-1]":500,
        "Separator heat conductivity [W.m-2.K-1]":0.1,
        "Positive particle radius [m]":10e-6,
        "Positive particle density [kg.m-3]":2500,
        "Positive particle heat capacity [J.kg-1.K-1]":700,
        "Positive particle heat conductivity [W.m-2.K-1]":2.1,
    }
)
r = pybamm.SpatialVariable(
    "r", 
    domain=["positive particle"],     
    auxiliary_domains={
        "secondary": "positive electrode"
    },
    coord_sys="spherical polar")
x_s = pybamm.SpatialVariable("x_s", domain=["separator"], coord_sys="cartesian")
x_p = pybamm.SpatialVariable("x_p", domain=["positive electrode"], coord_sys="cartesian")


geometry = {
    "separator": {x_s: {"min": -L_s, "max": 0}},        
    "positive electrode": {x_p: {"min": 0, "max": L_p}},
    "positive particle": {r: {"min": 0, "max": R_p}},
}
param.process_model(model)
param.process_geometry(geometry)

submesh_types = {
    "separator": pybamm.Uniform1DSubMesh,
    "positive electrode": pybamm.Uniform1DSubMesh,
    "positive particle": pybamm.Uniform1DSubMesh,
}
var_pts = {x_s: 10, x_p: 20, r: 30}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

spatial_methods = {
    "separator": pybamm.FiniteVolume(),
    "positive electrode": pybamm.FiniteVolume(),
    "positive particle": pybamm.FiniteVolume(),
}
disc = pybamm.Discretisation(mesh, spatial_methods)
disc.process_model(model);

[brosaplanella]

The issue is with the bc:

Ts: {"left": (pybamm.Scalar(0), "Neumann"), "right": (grad_Tb_cathode_particle_right, "Neumann")}

The problem is that the gradient evaluates on the edges of the discretised volumes (positive electrode in this case) but the particles are evaluated in the nodes of the positive electrode, so there's a mismatch. A possible solutions is to define grad_Tb_cathode_particle_right as

v = pybamm.PrimaryBroadcastToEdges(1, "positive electrode")
grad_Tb_cathode_particle_right = lumbda_n_eff / lumbda_s * pybamm.div(T_e_p * v)

which for the 1D case is the same you had before, but now it evaluates on the nodes. This is a bit of a hack, so probably there is a cleaner solution, but I am not familiar with it.

[xyyrly]

Thank you so much:)))This solution looks pretty good to me. I will try it~~