# Heat equation in non-dimensional form behaving differently than in usual format

Starting from $$c_p \frac{\partial u }{\partial t} = k \nabla^2 u$$ in a one dimensional domain [0,1] where $$c_p$$ and $$k$$ are modeling two different materials: $$k = \begin{cases} 1 ~\text{if} ~x < 0.5\\ 2.0 ~\text{else} \end{cases}$$ $$c_p = \begin{cases} 10^{-8} ~\text{if} ~x < 0.5\\ 1.0 ~\text{else} \end{cases}$$

I decided to refactor $$c_p$$ to the right hand side such that $$\frac{\partial u }{\partial t} = \frac{k} {c_p}\nabla^2 u$$ I solve both approaches using finite elements with Lagrange interpolation and Crank-Nicolson. However, both solutions are different (this is a random time step, the trend is similar for all time steps):

The one with $$c_p$$ refactored shows a flat solution for $$x<0.5$$, whereas the original equation has a linear solution. This difference disappears when the material properties are homogeneous, which makes me think I might be committing some mistake in my finite element formulation. The code to run both examples is:

from fenics import *

cp_electrolyte = 1e-8
k_electrolyte = 1.0
k_electrode = 2.0
cp_electrode = 1.0
scan_rate = 1.0
output_dir = "./"

mesh = UnitIntervalMesh(100)

V = FunctionSpace(mesh, "CG", 1)
u, v = TrialFunction(V), TestFunction(V)

Vlimit = 1.0
tlimit = Vlimit / abs(scan_rate)

class Materials(UserExpression):
def __init__(self, electrode, electrolyte, **kwargs):
super().__init__(**kwargs)  # This part is new!
self.electrolyte = electrolyte
self.electrode = electrode

def eval(self, values, x):
if x[0] < 0.5:
values[0] = self.electrolyte
else:
values[0] = self.electrode

k = Materials(k_electrode, k_electrolyte)
cp = Materials(cp_electrode, cp_electrolyte)

normal = False

def forward():

dt_value = 1e-2
dt = Constant(dt_value)
u_n = Function(V)
if normal:
a = cp * u / dt * v * dx + k * \
inner(Constant(1.0 / 2.0) * grad(u), grad(v)) * dx
L = (
cp * u_n / dt * v * dx
- k * inner(Constant(1.0 / 2.0) * grad(u_n), grad(v)) * dx
)
else:
a = u / dt * v * dx + k / cp * \
inner(Constant(1.0 / 2.0) * grad(u), grad(v)) * dx
L = (
u_n / dt * v * dx
- k / cp * inner(Constant(1.0 / 2.0) * grad(u_n), grad(v)) * dx
)

t = 0
T = tlimit * 5
n_steps = int(T / dt_value)

bcval = Expression("t", t=t, degree=1)

def Left(x, on_boundary):
return x[0] < DOLFIN_EPS and on_boundary
bc = DirichletBC(V, bcval, Left)

u_sol = Function(V)
if normal:
output = "potential.pvd"
else:
output = "potential_ratio.pvd"
potential_pvd = File(output)
while t < T:
solve(a == L, u_sol, bcs=bc)
t += dt_value
bcval.t = t
potential_pvd << u_sol
u_n.assign(u_sol)

return u_n

u_n = forward()



Thanks

## 2 Answers

The applied division is fine, what went wrong here, is the application of Stoke's theorem. If you multiply with the test function you get following term:

$$\int \frac{1}{c_p}\nabla\left(-k\nabla u\right) v d\Omega$$ But
$$\int \frac{1}{c_p}\nabla\left(-k\nabla u\right) v d\Omega \neq \int \frac{1}{c_p} \left(k\nabla u\right) \cdot \left(\nabla v\right) d\Omega+\int \frac{1}{c_p} (-k \nabla u) v d\partial\Omega$$ You miss the derivative of $$\frac{1}{c_p}$$, exactly: $$-\int \left(\nabla \frac{1}{c_p}\right) \cdot \left(-k\nabla u\right) v d\Omega$$

If you add this term in the second formulation, the results should be identical.

One more thing, the heat flux $$q=-k\nabla u$$, so if $$k$$ is dependent on $$x$$, you should not move it outside of the divergence operator. Also for completeness you miss $$\rho$$ in front of the time derivative as well.

Edit

While I agree that integrating $$c_p$$ inside the mass matrix is much, much simpler, the reason why the second variant fails, is the missing derivative of the function $$\frac{1}{c_p}$$. This missing description of the change in value in $$c_p$$ is present if used in the mass matrix.

So how can you incorporate this?

$$c_p$$ can be written as $$c_p=a + H(x-x_0) (b-a)$$ with $$H(\cdot)$$ the Heavyside step function which has as derivative $$\delta(x)$$ the Dirac delta function. $$a, b$$ and $$x_0$$ the corresponding values.

Now, the contribution $$\nabla(\frac{1}{c_p})=\frac{-1}{c_p^2} (b-a) \delta(x-x_0)$$ is $$-\int \left(\frac{-1}{c_p^2} (b-a) \delta(x-x_0)\right) \cdot \left(-k\nabla u\right) v d\Omega \neq 0$$.

The simplest way to incorporate such a contribution in finite element methods is to use a mollified version of the distribution. I am not an expert on fenics, but here is an example of this. It is likely a good idea to use a mollifier on the Heavyside step function and calculate its derivative to make it consistent.

Second Edit

The following code demonstrates that the error is indeed in the application of Stoke's theorem.

import fenics
import gmsh
import numpy
from fenics import (DOLFIN_EPS, Constant, DirichletBC, Expression, File,
Function, FunctionSpace, TestFunction, TrialFunction,
UnitIntervalMesh, UserExpression, atan, dx, grad, inner,
interpolate, nabla_grad, pi, project, solve)

cp_electrolyte = 1e-8
k_electrolyte = 1.0
k_electrode = 2.0
cp_electrode = 1.0
scan_rate = 1.0
output_dir = "./"
eps=1e-4

mesh = UnitIntervalMesh(1000)

V = FunctionSpace(mesh, "CG", 1)
Q = FunctionSpace(mesh, "DG", 0)  # this space for the derivative of cp
u, v = TrialFunction(V), TestFunction(V)

Vlimit = 1.0
tlimit = Vlimit / abs(scan_rate)

def heavyside(x, eps):
value = 1 / 2 + 1 / pi * atan(x / eps)
return value

class Materials(UserExpression):
def __init__(self, electrode, electrolyte, **kwargs):
super().__init__(**kwargs)  # This part is new!
self.electrolyte = electrolyte
self.electrode = electrode

def eval(self, values, x):
if x[0] < 0.5:
values[0] = self.electrolyte
else:
values[0] = self.electrode

class MaterialsEps(UserExpression):
def __init__(self, electrode, electrolyte, eps, **kwargs):
super().__init__(**kwargs)
self.electrolyte = electrolyte
self.electrode = electrode
self.eps = eps

def eval(self, values, x):
values[0] = self.electrolyte + \
heavyside(x[0]-0.5, self.eps)*(self.electrode-self.electrolyte)

k = Materials(k_electrode, k_electrolyte)

# show material functions
cp = project(Materials(cp_electrode, cp_electrolyte), V)
mat = File("material.pvd")
cp.rename("projected material function", "")
mat << cp

# mollified material function
cp = project(MaterialsEps(cp_electrode, cp_electrolyte, eps), V)
dcpdx = project(cp.dx(0), Q)

mat = File("materialeps.pvd")
cp.rename("projected mollified material function", "")
mat << cp

mat = File("gradmaterialeps.pvd")
dcpdx.rename("projected derivative of mollified material function", "")
mat << dcpdx

def forward(eps):
dt_value = 1e-2
dt = Constant(dt_value)
u_n = Function(V)
heatflux = -1 * k * grad(u)
heatflux_n = -1 * k * grad(u_n)

if normal:
cp = Materials(
cp_electrode,
cp_electrolyte)  # original discontinous material function
a = cp * u / dt * v * dx - 1 / 2 * (inner(heatflux, grad(v)) * dx)
L = cp * u_n / dt * v * dx + 1 / 2 * (inner(heatflux_n, grad(v)) * dx)
else:
# mollified material function
cp = project(MaterialsEps(cp_electrode, cp_electrolyte, eps), V)
dcpdx = project(cp.dx(0), Q)
a = u * v * dx - 1 / 2 * dt * (
inner(heatflux, grad(v)) / cp * dx + 1 /
(cp * cp) * dcpdx * k * nabla_grad(u)[0] * v * dx)
# with missing derivative term
# likely, there is a better solution than explicitely access the first element of nabla_grad
L = u_n * v * dx + 1 / 2 * dt * (
inner(heatflux_n, grad(v)) / cp * dx + 1 /
(cp * cp) * dcpdx * k * nabla_grad(u_n)[0] * v * dx)

t = 0
T = tlimit * 5
n_steps = int(T / dt_value)

bcval = Expression("t", t=t, degree=1)

def Left(x, on_boundary):
return x[0] < DOLFIN_EPS and on_boundary

bc = DirichletBC(V, bcval, Left)

u_sol = Function(V)
if normal:
output = "potential.pvd"
u_sol.rename("u_original", "")
else:
output = "potential_ratio.pvd"
u_sol.rename("u_with_mollified_term", "")
potential_pvd = File(output)
while t < T:
solve(a == L, u_sol, bcs=bc)
t += dt_value
bcval.t = t
potential_pvd << u_sol
u_n.assign(u_sol)

return u_n

for cc in [True, False]:
normal = cc
u_n = forward(eps)


• $c_p$ is a function of $x$ and is discontinuous. How do you propose to do integration by parts in this case ? Formally you get a volume term like $\int k \nabla u \cdot \nabla (v/c_p) dx$. How do you propose to compute derivatives involving $c_p$ ? – cfdlab Oct 22 '20 at 3:20
• @cfdlab see edit. – Bort Oct 22 '20 at 9:41
• Thanks! I wonder if the integral with the dirac delta could be represented with a line integral (point integral in the 1D case). Although I do not know what value of $1 / c_p^2$ at that point one should take. – balborian Oct 22 '20 at 15:56
• @balborian does my posted code answers your question or are you wondering how to do this with a direct manipulation of the matrix? – Bort Oct 22 '20 at 16:12
• It does answer my question, but I am wondering if there is a cleaner way to do it rather than using the approximation of the dirac delta. – balborian Oct 22 '20 at 16:14

If $$k$$ depends on the spatial variables, the heat equation is of the form $$c_p u_t = \nabla \cdot (k \nabla u)$$ In your case, $$c_p$$ also depends on space and is discontinuous. You should not try to divide by $$c_p$$ in this case. At best you can divide by some constant $$c_{p,ref}$$, e.g. $$c_{p,ref} = \max_{x,y} c_p(x,y)$$ which is just one in your case. You have to integrate $$c_p$$ inside the mass matrix. The semi-discrete scheme would look like this $$M du/dt = A u$$ $$M_{ij} = \int c_p \phi_i \phi_j dx, \qquad A_{ij} = -\int k \nabla\phi_i \cdot \nabla \phi_j dx$$ together with some boundary conditions. You can put in a time discretization after this.