How to solve transient equations?

[CalebDmArcher]

Hi, I need to solve the heat diffusion equation shown below:

$$ \nabla \cdot \left( k(x, y, z) \nabla T(x, y, z) \right) + \dot{q}(x, y, z) = \rho c_p \frac{\partial T}{\partial t} $$

The equation I am having right now is:

eq_transient = TransientTerm(coeff=rho*cp) == (DiffusionTerm(coeff=Gamma) +
                                               source +
                                               (RobinCoeff * g).divergence -
                                               ImplicitSourceTerm(coeff=(RobinCoeff * a.dot(n)).divergence))

where rho is the density and cp is the absolute heat cap in the heat diffusion equation.
I am not sure if this is the correct way to write the transient term.

The way I am using to solve the equation is shown below:

timeStep = 1.e-7
steps = 50
for step in range(steps):
    eq_transient.solve(var=temperature, dt=timeStep)
    print('Current step number is: ', step)

I noticed that as I decreased the timeStep from 1.e-1 to 1.e-7 with 'steps' values from 50 to 300, the final results converged to a certain value shown in the table below. (for the same timeStep value, the results are the same with different steps values)

I was doing steady-state equation solving, but the steady-state results are very weird. Such as when there is no heat source, part of the material has a temperature of 0 °C while the ambient temperature is 45 °C; when there is a heat source, the temperature can go over 10000 °C (while the transient result is around 45 °C shown in the table above). So I think I might need to try to do transient analysis with proper time steps to reach the true steady state results.

My questions are:

1. Is this the correct way of writing the transient equation?
2. Is this the correct way to solve the transient equation? If so, do the different results from different timeStep mean that I need to figure out the best timeStep for different simulation setups?
3. Is the unit of timeStep second OR how can I know the unit / exact length of it?

The full code is attached below:
example01_07_25.zip

[guyer]

1. Is this the correct way of writing the transient equation?

Sure.

2. Is this the correct way to solve the transient equation? If so, do the different results from different timeStep mean that I need to figure out the best timeStep for different simulation setups?

A transient problem needs to call temperature.updateOld() in order to advance the time step. As it is, your loop is just solving the same time step over and over again. Different time steps give different answers because they have evolved for different amounts of time, mostly nowhere near "steady state" [*].

3. Is the unit of timeStep second OR how can I know the unit / exact length of it?

The unit of timeStep is whatever is consistent with the units of your coefficients. You give $\rho$ as [g/m^3], $c_p$ as [J/(g*K)], $k$ as [W/(m*K)], and the dimensions of your mesh in meters, so the units of timeStep are [J/W], i.e., [s].

As before, I strongly recommend you simplify your problem and solve it in lower dimensions. This problem is fundamentally 1D; if you change to a Grid1D and add a Viewer, it's much easier to see what's going on.

[*] If you do that, what you can see is that you don't really ever get a proper steady-state, although it does seem to evolve to a steady shape with a steady $\Delta T$. At least one reason is that you have a constant heat source pumping into the system, so the temperature rises without bound. A further issue is that this equation admits an infinite number of solutions. Without any boundary condition to "ground" the solution, it may tend to float, although the Robin condition should keep things near T_ambient (if not for the source).

[CalebDmArcher]

Thank you so much for the reply. Hope you had a great Christmas and New Year.

I will try to write the code in 1D. Will let you know if I can get something.

I did a calculation by hand with only the free air convection on the top and the bottom surface and tried to solve the top surface temperature by using heat flow and thermal resistors. An approximation is made, which assumes all power is applied to the top surface. And intuitively thing, I feel the results from the heat-on-surface case should be higher than the volume-heat-generation case since now power is concentrated on the surface so I guess the surface should be warmer. The derivation is shown below
(There is a mistake in this image so I deleted it.)
and the results with different q values are:
(There is a mistake in this image so I deleted it.)
Clearly, there is a discrepancy between what I thought and what the 3D simulation results showed (after adding temperature.updateOld() to the time step for loop, I did see a temperature increasing way above 45 as the steps value increasing and also consistency with period, such as the maximum value of timeStep = 1.e-8 running steps = 500 gets the same value as timeStep = 1.e-7 running steps = 50).

I will work on 1D and see if I can spot any mistakes in my code.
Sorry for being annoying. Really appreciate your help and your help really means a lot to me🥹.

[CalebDmArcher]

Update:
I spotted a mistake in the code: I missed the dot in this line k_si = 149. so the system had weird results previously.
The theoretical value calculated by hand is

and with different heat source values:

After fixing this mistake, the steady-state result matches the hand-calculate result.