Planar transient heat problem

Compute the transient temperature distribution in a planar problem.

import matplotlib.pyplot as plt
import torch
from matplotlib import animation

from torchfem.data import get_data
from torchfem.io import import_mesh
from torchfem.materials import IsotropicConductivity2D
from torchfem.planar import PlanarHeat
from torchfem.plot_utils import embed_animation_gif

# Set default data type to double precision
torch.set_default_dtype(torch.float64)

# Dimensions (thickness of plate t and edge length L)
t = 0.01
L = 0.2

Transient heat equation on planar domain

# material properties
kappa = 500.0
rho = 7850.0
cp = 500.0

# Material model (thermal mass is rho*cp)
material = IsotropicConductivity2D(kappa=kappa, rho=rho * cp)

# Import mesh
mesh = import_mesh(get_data("plate_hole.vtk"), material)

# origin in center
mesh.nodes[:, :2] -= 0.5 * L

# get radius
distance_to_center = torch.linalg.norm(mesh.nodes, dim=-1)
R = distance_to_center.min()

plate = PlanarHeat(mesh.nodes, mesh.elements, material)
plate.thickness[:] = t

# Fixed temperature at left end
left = plate.nodes[:, 0] < -0.5 * L + 1e-6
plate.constraints[left] = True
plate.displacements[left, 0] = 5.0

# Fixed temperature at right end
right = plate.nodes[:, 0] > 0.5 * L - 1e-6
plate.constraints[right] = True
plate.displacements[right, 0] = 20.0

# Fixed temperature at radius
radius = torch.isclose(distance_to_center, R)
plate.constraints[radius] = True
plate.displacements[radius, 0] = 0.0


# solver settings
end_time = 40.0
delta_t = 0.1

# Results are returned at exactly these times
times = torch.arange(0.0, end_time + delta_t, delta_t)

# Solve
temp, rfl, hf, temp_grad, _ = plate.time_integration(times, delta_t, verbose=False)

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(10, 10))


def plot_state(i):
    plate.plot(
        node_property=temp[i],
        cmap="magma",
        title="Temperature",
        ax=ax1,
        vmin=0.0,
        vmax=20.0,
    )
    plate.plot(
        node_property=rfl[i],
        cmap="coolwarm",
        title="Surface heat flux",
        ax=ax2,
    )
    plate.plot(
        element_property=-hf[i],  # negative sign
        cmap="coolwarm",
        title="Heat flux",
        color="lightgray",
        ax=ax3,
    )
    plate.plot(
        element_property=temp_grad[i],
        title="Temperature Gradient",
        ax=ax4,
        cmap="viridis",
        color="lightgray",
    )

    # if i == 0:
    #     plt.tight_layout()


def animate(i):
    # Note. this is pretty slow and should be optimized.
    # updating the data on the tricontour is not intended in mpl.
    ax1.clear()
    ax2.clear()
    ax3.clear()
    ax4.clear()

    plot_state(i)

    return ()  # we have to return a (empty) tuple of artists


anim = animation.FuncAnimation(
    fig,
    animate,
    init_func=plot_state(0),
    frames=range(len(temp))[::10],
    interval=20,
    blit=False,
    repeat=False,
)

plt.close(fig)  # prevent static plot
embed_animation_gif(anim, fps=10)  # <- This embeds the GIF directly