Problem simulating heat transport through composite solid

Hello NGSolve users,

I have a problem simulating the heat transport through my model of a stirling engine, which is a solid composed of different materials.
I modeled the geometry in FreeCad, exported the .step file, created a mesh using GMSH and used the mesh in the NGSolve script using ReadGmsh().
I added the .step, .geo and .msh files as well as the simulation script in the attachements.

In my Simulation, the lower surface is set to a constant temperature higher than the initial temperature, such that the body is heating up from below. The problem is, that there are some artefacts in the temperature field, and i struggle to find out whether the mesh file is corrupted or the python script is faulty (I think i already checked the second option using another geometry and mesh file). The temperature looks like this
and thus illogical to me.
I would be very grateful if someone could take a look at my code and give me a hint.

Kind regards and thanks in advance
Lukas

Hi Luma,

when you have different materials, you have to keep the densities and heat capacities with the mass matrix, otherwise the interface conditions at the material boundaries are wrong.

I am not sure how you got to your implicit Euler time-stepping. I replaced it by my standard one.

The attached file gives reasonable results.

Best, Joachim

Hi Joachim,

thank you for the quick reply, I have been unaware of that. Works perfectly fine now!

Have a nice day
Lukas

I must correct:
If i set the timestep and end time back to dt=0.01 and t_end=1, the result is the same as before. You have any idea why that is?

And as for keeping the densities and heat capacities with the mass matrix, i dont really understand why it makes a difference:
Solving
[tex]
\frac{\partial}{\partial t}u - \frac{\lambda}{\rho c}\Delta u = f
[/tex]
and
[tex]
\frac{\partial}{\partial t}\rho c u - \lambda\Delta u = \rho c f
[/tex]
should yield the same result if i set the diffusivity right on the different subdomains.

now you see an underresolved initial layer which occurs since initial conditions and boundary conditions don't fit together. just wait.

with the weak form we discretize
div (lambda grad u)

if lambda is not one constant everywhere, you cannot move it through the div operator.