Scikit-fem wire heat simulation problem (Crank-Nicolson)

[Noekil]

I'm trying to simulate the temperature evolution in a cable due to the Joule effect. I already have a code that allows me to simulate conduction in a cable but not temporally. I've heard that one possible method is the Crank-Nicolson method. I've tried to look into it but I can't seem to understand how to apply this method to my code. Maybe there's an expert in this kind of simulation who will come across this message. Thanks in advance for the help!

from pathlib import Path
from typing import Optional
import numpy as np

from skfem import *
from skfem.helpers import dot, grad
from skfem.models.poisson import mass, unit_load
from skfem.io.json import from_file

joule_heating = 5.
heat_transfer_coefficient = 7.
thermal_conductivity = {'core': 101., 'annulus': 11.}

mesh = from_file(Path(file).parent / 'meshes' / 'disk.json')
radii = sorted([np.linalg.norm(mesh.p[:, mesh.t[:, s]], axis=0).max() for s in mesh.subdomains.values()])

@BilinearForm
def conduction(u, v, w):
return dot(w['conductivity'] * grad(u), grad(v))

convection = mass

element = ElementQuad1()
basis = Basis(mesh, element)

basis0 = basis.with_element(ElementQuad0())

conductivity = basis0.zeros()
for s in mesh.subdomains:
conductivity[basis0.get_dofs(elements=s)] = thermal_conductivity[s]

L = asm(conduction, basis, conductivity=basis0.interpolate(conductivity))

facet_basis = FacetBasis(mesh, element, facets=mesh.boundaries['perimeter'])
H = heat_transfer_coefficient * asm(convection, facet_basis)

core_basis = Basis(mesh, basis.elem, elements=mesh.subdomains['core'])
f = joule_heating * asm(unit_load, core_basis)

temperature = solve(L + H, f)

T0 = {
"skfem": (basis.probes(np.zeros((2, 1))) @ temperature)[0],
"exact": (
joule_heating
* radii[0] ** 2
/ 4
/ thermal_conductivity["core"]
* (
2 * thermal_conductivity["core"] / radii[1] / heat_transfer_coefficient
+ (
2
* thermal_conductivity["core"]
/ thermal_conductivity["annulus"]
* np.log(radii[1] / radii[0])
)
+ 1
)
)
}

if name == 'main':

from os.path import splitext
from sys import argv
from skfem.visuals.matplotlib import draw, plot, savefig

ax = draw(mesh)
plot(basis, temperature, ax=ax, colorbar=True)
savefig(splitext(argv[0])[0] + '_solution.png')

I hope someone can explain to me how to make this simulation temporal. If anyone has an example, please don't hesitate to share.

[kinnala]

If you discretize the heat equation

$u' = \Delta u$

in space using FEM, you will get the system of ODEs:

$M x' = Kx.$

Then you can apply finite difference method in time using time step $\delta$, e.g., implicit Euler method:

$M \frac{x_k - x_{k-1}}{\delta} = K x_k.$

You can rearrange the equation to get a linear system for $x_k$:

$(M - \delta K) x_k = M x_{k-1}.$

Then in code you implement matrices $M$ and $K$, and do a for-loop over time steps with an initial guess for $x_0$.

[gdmcbain]

Here are examples of solving the heat equation using Crank–Nicolson:

• ex19: heat equation
• ex39: one-dimensional heat equation

These use a generalization of the Crank–Nicolson method called the θ-method which reduces to Crank–Nicolson for θ = ½. It's derived pretty much as by @kinnala (above) for the special case θ = 1; It also gives the explicit Euler method for θ = 0.