The Weak Form of Heat Conduction

This section describes the physics of steady-state heat conduction used in Scikit-Topt. We derive the governing equation, weak formulation, and treatment of Dirichlet, Neumann, and Robin boundary conditions in the finite element method (FEM).

Governing Equation (Strong Form)

Let \(T : \Omega \rightarrow \mathbb{R}\) be the temperature field. The steady-state heat conduction equation without internal heat source is:

\[-\nabla \cdot (k(\mathbf{x}) \nabla T) = 0 \quad \text{in } \Omega,\]

where

  • \(k(\mathbf{x})\) is the thermal conductivity

  • \(T\) is the unknown temperature

Here, the thermal conductivity \(k(\mathbf{x})\) is treated as a scalar (isotropic conductivity). Anisotropic or tensor-valued conductivity is not considered in the current implementation.

Boundary Conditions

We consider three types of thermal boundary conditions:

  1. Dirichlet (Temperature BC) Prescribed temperature:

    \[T = T_0 \quad \text{on } \Gamma_D.\]

  2. Neumann (Heat Flux BC) Prescribed heat flux \(q\):

    \[-k \nabla T \cdot \mathbf{n} = q \quad \text{on } \Gamma_N.\]

  3. Robin (Convective BC) Newton’s cooling law:

    \[-k \nabla T \cdot \mathbf{n} = h(T - T_\infty) \quad \text{on } \Gamma_R,\]

    where

    • \(h\) : heat transfer coefficient

    • \(T_\infty\) : ambient temperature

    This boundary condition models cooling by air or fluid.

Weak Formulation

Let \(v\) be a test function. Multiply the strong form by \(v\) and integrate over \(\Omega\):

\[\int_\Omega -\nabla \cdot (k \nabla T) \, v \, d\Omega = 0.\]

Using the divergence theorem:

\[\int_\Omega k \nabla T \cdot \nabla v \, d\Omega = \int_{\partial\Omega} (k \nabla T \cdot \mathbf{n}) v \, d\Gamma.\]

We now insert boundary conditions.

Dirichlet BC

For \(\Gamma_D\), the temperature is prescribed. As usual in FEM, this is enforced through modification of the linear system after assembling the matrix (see below).

Neumann BC

On \(\Gamma_N\):

\[-k \nabla T \cdot \mathbf{n} = q.\]

Substituting into the boundary integral:

\[\int_{\Gamma_N} q v \, d\Gamma.\]

Robin BC

Robin BC provides an additional contribution:

\[-k \nabla T \cdot \mathbf{n} = h(T - T_\infty) \quad \Rightarrow \quad k \nabla T \cdot \mathbf{n} = -h(T - T_\infty).\]

Plug into the boundary integral:

\[\int_{\Gamma_R} h(T - T_\infty) v \, d\Gamma.\]

This term splits into:

\[\int_{\Gamma_R} h T v \, d\Gamma - \int_{\Gamma_R} h T_\infty v \, d\Gamma.\]

As a result:

  • The first part contributes to the matrix (Robin = “boundary stiffness”)

    \[K_{ij}^{(R)} = \int_{\Gamma_R} h \phi_i \phi_j \, d\Gamma.\]

  • The second part contributes to the load vector

    \[f_i^{(R)} = \int_{\Gamma_R} h T_\infty \phi_i \, d\Gamma.\]

FEM Discretization

Expand the temperature field as:

\[T_h(\mathbf{x}) = \sum_{j=1}^N T_j \phi_j(\mathbf{x}).\]

The weak form becomes:

\[\mathbf{K}\mathbf{T} = \mathbf{f},\]

where

  • stiffness matrix

    \[K_{ij} = \int_\Omega k \nabla\phi_j \cdot \nabla\phi_i \, d\Omega + \int_{\Gamma_R} h \phi_j \phi_i \, d\Gamma.\]

  • load vector

    \[f_i = \int_{\Gamma_N} q \phi_i \, d\Gamma + \int_{\Gamma_R} h T_\infty \phi_i \, d\Gamma.\]

Dirichlet Boundary Enforcement

Dirichlet temperature BCs (prescribed \(T=T_0\)) are applied by modifying the linear system:

  • Zero out rows of fixed DOFs

  • Set diag = 1

  • Set load entry = prescribed value

In Scikit-FEM (as used by Scikit-Topt):

Solving the Linear System

The final assembled system is:

\[\mathbf{K}_e \mathbf{T} = \mathbf{f}_e.\]

In Scikit-Topt, the solver is typically:

T = scipy.sparse.linalg.spsolve(K_e, f_e)