About
Scikit-fem is a lightweight and pure-Python finite element library designed for flexibility, transparency, and educational clarity. Unlike large compiled frameworks, it allows users to express weak forms directly in Python syntax, making it easy to prototype new methods or couple physics without the complexity of external dependencies. Despite its simplicity, Scikit-fem supports multi-dimensional meshes, custom element definitions, and efficient sparse-matrix assembly through NumPy and SciPy backends. Its modular design makes it particularly suitable for research, teaching, and the rapid development of custom finite element solvers.
New Feature
Building upon Scikit-fem, I implemented an additional module to solve steady-state heat conduction problems governed by the equation
\nabla \cdot (k \nabla T) = q
with Dirichlet, Neumann, and Robin boundary conditions.
This formulation enables the simulation of heat transfer in solid domains while accounting for convection and prescribed temperature boundaries.
Furthermore, I extended this framework to perform thermal topology optimization by minimizing thermal compliance, defined as
C = \int_\Omega q T d\Omega
subject to a volume constraint on the material distribution.
Using a SIMP(Solid Isotropic Material with Penalty)-type interpolation for the thermal conductivity
k(\rho) = k_{\min} + \rho^p (k_0 - k_{\min})
the optimization loop updates the density field ρ\rhoρ to design structures that efficiently conduct heat away from hot regions.
This extension demonstrates how Scikit-fem’s flexible weak-form interface can be leveraged to implement advanced multi-physics optimization frameworks entirely in Python.
And lastly, the drawback of using the SIMP method and possible directions for future work are discussed.
Topology Optimization
Task Definition
The topology optimization problem was defined by setting the green region as the design domain, the blue region as the heat source (Dirichlet boundary condition), and the red region as the Robin boundary condition.
No other boundary conditions were applied, so heat dissipation occurs only through the Robin boundaries.
The Robin boundary condition is an intermediate condition between the Dirichlet (prescribed temperature) and Neumann (prescribed heat flux) boundary conditions, and is commonly used to model heat transfer and thermal dissipation.
-k \frac{\partial T}{\partial n} = h (T - T_{\mathrm{env}}) \quad \text{on } \Gamma_RHere,
- ( \displaystyle \frac{\partial T}{\partial n} = \nabla T \cdot n ): temperature gradient in the outward normal direction
- ( h ): heat transfer coefficient (convection coefficient)
- ( T_{\mathrm{env}} ): ambient or environmental temperature
- ( \Gamma_R ): Robin boundary surface
You can also see below article here,
The thermal compliance was defined as described above, and optimization was performed using the SIMP-based density gradient method.
The compliance is extended to
\int_\Omega q T d\Omega = \int_\Omega k(\rho) |\nabla T|^2 d\Omega.
To derive the derivatives from the compliance,
\frac{\partial C}{\partial \rho_e}
= - \int_{\Omega_e} \frac{\partial k(\rho_e)}{\partial \rho_e} |\nabla T|^2 d\Omega.SIMP interpolation is used above, it is described as
k(\rho_e) = k_{\min} + \rho_e^p (k_0 - k_{\min})therefore,
\frac{\partial k(\rho_e)}{\partial \rho_e}
= p \rho_e^{p-1} (k_0 - k_{\min})\\
\frac{\partial C}{\partial \rho_e}
= - p \rho_e^{p-1} (k_0 - k_{\min})
\int_{\Omega_e} |\nabla T|^2 , d\Omega.
Result
The drawback of SIMP and Future Work
What I realized after performing thermal topology optimization using the SIMP method is that it is ultimately difficult to obtain an ideal boundary configuration with SIMP.
This is because the SIMP method performs optimization based on predefined Robin boundary conditions—in other words, the heat-dissipating surfaces must be specified in advance.
As a result, heat dissipation along intermediate internal surfaces that naturally emerge during the optimization process is not taken into account.
There are several studies attempting to address this issue, but in the end, level-set–based optimization using unfitted meshes may provide a more appropriate framework—although its implementation is considerably more complex.
Let’s modify the previous problem slightly as a test.
In the earlier example, three separate regions at the top were assigned Robin boundary conditions.
This time, we will assign the Robin condition to the entire top surface instead.
As a result, we obtain a rather uninteresting density distribution, as shown below.
Why does this happen?
It is because the Robin boundary condition is not directly applied to the boundaries generated by the SIMP method.
As a result, no heat dissipation occurs along those newly formed interfaces.
Consequently, the optimizer tends to create straight conductive paths that connect the heat source directly to the predefined Robin boundary.
In structural optimization, this is not a serious issue, since there is no need to impose penalties on boundaries that emerge during the optimization process.
However, in thermal problems, introducing boundary-related penalties or heat-loss mechanisms is advantageous—something that the SIMP framework is not naturally designed to handle.
Of course, there have been several studies attempting to overcome this limitation.
For example:
- Introducing penalty terms similar to the Robin boundary condition on internal interfaces,
- Making the heat transfer coefficient of the Robin condition density-dependent, or
- Employing alternative loss formulations to better capture heat dissipation effects.
Insight from Other works
Since the heat conduction equation is almost identical to that of structural analysis, I implemented a nearly equivalent compliance optimization. However, according to other papers, it appears that the objective function for heat sink optimization is the average temperature.
See material below:
The objective of them is different though, they use the average temprature.
So what is the difference from compliance?
The objective of compliance
Its objective is to make the overall heat flux more uniform and reduce gradients.
→ It tends to smooth out the temperature distribution.
However, it does not necessarily create the shortest or most direct path to cool the heat source.
The objective of average temperature
Its objective is to cool the heat-generating region.
→ It tends to create efficient paths around the heat source to dissipate heat.
Future Work
It’s interesting that minimizing the average temperature is used with the purpose of lowering the temperature of the heat source itself.
Considering applications such as heat sink design, this approach indeed seems more appropriate.
Then, what kinds of industrial applications might benefit from heat compliance minimization?
Since reducing temperature nonuniformity helps prevent thermal stress, deformation, and material degradation, it could be useful for mold, casting, and injection molding processes where suppressing thermal stress, distortion, and fatigue is important.
It might also be applicable to circuit boards and similar components.
In any case, the next task will likely be to implement the objective function and derivative for average temperature minimization.
