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

However, both solutions with FEniCS are different (this is a random time step, the trend is similar for all time steps):
image

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

Note that the weak form \int_\Omega k\nabla u\cdot\nabla v\,d\mathbf{x} corresponds to the strong form -\nabla\cdot(k\nabla u), which is not the same as -k\nabla^2 u when k varies in space. Likewise, when c_p^{-1} varies in space, it cannot be pulled inside the divergence operator, so

-c_p^{-1}\nabla\cdot(k\nabla u) \neq -\nabla\cdot(c_p^{-1}k\nabla u)\text{ ,}

where the left side is what you get by dividing the original conservation law through by c, but the right side is what the weak form in the code corresponds to when normal == False.

While k and c_p^{-1} are piecewise-constant, there are Dirac deltas in their spatial derivatives at the discontinuity. (You could pull c_p^{-1} inside the divergence operator on either side of the internal interface where c_p and k are discontinuous, but then when you integrate by parts on either side to get a weak form, there would be non-cancelling boundary terms at the interface.)

1 Like

For anyone interested, someone posted the code in scicomp StackExchange: https://scicomp.stackexchange.com/a/36143/7023 . The only issue is that they use an approximation to the dirac delta.