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

Question

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

Answer 1

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)

Answer 2

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.