I am trying to solve a non-linear time-dependent heat equation
$$\partial_tT=\nabla \left(k_T(T)\nabla T\right) + f$$ (similar to question Solving a non-linear heat equation with the galerkin method gives negative values) with a FEM-approach, using the continuous Galerkin method for the FEM part, the Newton method for linearization of the non-linear part and the implicit Euler method for the time stepping.
Initially, I had problems with parts of the results becoming negative, which could be solved (more or less) by clipping all parts which are below a certain value. Now I got the additional problem that the solution at some nodes becomes unstable, and due to a positive feedback loop increases into unrealistic values. This leads to results as shown in the following image:
The heat source is located at the bottom, approximately at 500, which can be seen due to the darker yellow color. Still, in the middle of the image, the solution shows some nodes with values significantly higher values than the surrounding area, leading to a non-smooth result.
Thus I was wondering if there is a method to prevent those runaways, such as a slope limiter (and if yes, how large will the error be)? I can not use a lowpass here, after those values can be reached during regular calculations.
Or is my approach in general flawed for such an equation, and I should switch to another method?
To address the comment: In general I try to solve the equations from https://research.utwente.nl/en/publications/two-temperature-model-for-pulsed-laser-induced-subsurface-modific, with the equation for the electron temperature ($T_E$) resulting in node spikes, while the other two equations do not make any problems. The boundary conditions are Neumann conditions with $dT/dt=0$. The domain itself is defined as a square from (0, 0) to (1, 1), due to the problem being dimensionless, while the basis functions are polynomials.