Planar stationary heat problem

Compute the stationary temperature distribution in a planar plate with prescribed temperatures at its left and right edges and convection to ambient at its top edge.

import matplotlib.pyplot as plt
import torch

from torchfem.data import get_data
from torchfem.io import import_mesh
from torchfem.materials import IsotropicConductivity2D

# 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

Static heat equation on planar domain

# Material model
material = IsotropicConductivity2D(kappa=1000.0)

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

# Fixed boundary at left end
left = plate.nodes[:, 0] < 1e-6
plate.constraints[left] = True
plate.temperatures[left, 0] = 10.0

# Fixed boundary at right end
right = plate.nodes[:, 0] > L - 1e-6
plate.constraints[right] = True
plate.temperatures[right, 0] = 20.0

# Convection to ambient at top edge
top = plate.nodes[:, 1] > L - 1e-6
plate.convection(top, h=5000.0, T_inf=0.0)

# Solve
temp, rfl, hf, temp_grad, _ = plate.solve(verbose=True)

────────────────────────────────────────────────────────────────────────────────────────
 model    PlanarHeat · 790 elem · 451 dof · float64
 machine  x86_64 · 2 threads · 16 GB RAM
 solver   direct · scipy · cpu
 newton   rtol 1e-08 · atol 1e-06 · <=10 it
────────────────────────────────────────────────────────────────────────────────────────
  Increment  Load factor  Steps  Newton  Krylov  Residual  Wall time
          1            1      1       1       -  1.66e-12     0.01 s
────────────────────────────────────────────────────────────────────────────────────────
 converged · 1 increment · 1 Newton iteration · 0.01 s

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

plate.plot(
    node_property=temp,
    cmap="magma",
    title="Temperature",
    ax=ax1,
)
plate.plot(
    node_property=rfl,
    cmap="coolwarm",
    vmin=-rfl.abs().max(),
    vmax=rfl.abs().max(),
    title="Surface heat flux",
    ax=ax2,
)
plate.plot(
    element_property=hf,
    cmap="coolwarm",
    title="Heat flux",
    color="white",
    edgecolor="lightgray",
    ax=ax3,
)
plate.plot(
    element_property=temp_grad,
    title="Temperature Gradient",
    ax=ax4,
    cmap="viridis",
    color="white",
    edgecolor="lightgray",
)
plt.tight_layout()
plt.show()