Transient heat diffusion

The next step is to derive and implement the unsteady (time-dependent) heat diffusion equation. The transient heat equation looks like this:

(140)\[\rho c_p \frac{\partial T}{\partial t} = \frac{\partial}{\partial x} k \frac{\partial T}{\partial x} + \frac{\partial}{\partial y}k\frac{\partial T}{\partial y}\]

Two differences to the previously considered steady-state heat diffusion are apparent. First, we now have a time derivative in addition to the spatial derivatives. Second, the material parameters (density, specific heat, thermal conductivity) are no longer assumed constant. Nevertheless, the diffusion part of equation (140) looks very similar to the steady-state case that we know how to solve.

We know how to handle spatial derivatives but this is the first time we encounter a time derivative in finite elements. We use a first-order implicit (backward Euler) finite-difference for time:

(141)\[\rho c_p \frac{T^{n+1} - T^n}{\Delta t} = \frac{\partial}{\partial x} k\frac{\partial T^{n+1}}{\partial x} + \frac{\partial}{\partial y}k\frac{\partial T^{n+1}}{\partial y}.\]

Re-arrange equation (141) so that all known temperatures \(T^n\) are on the Rhs and all unknown temperatures \(T^{n+1}\) are on the Lhs.

(142)\[\rho c_p T^{n+1} - \Delta t \left( \frac{\partial}{\partial x} k\frac{\partial T^{n+1}}{\partial x} + \frac{\partial}{\partial y}k\frac{\partial T^{n+1}}{\partial y} \right ) = \rho c_p T^{n}.\]

The backward Euler scheme is unconditionally stable for linear diffusion (useful for larger time steps), but only first-order accurate in time.

FEM form

Now we proceed in the usual way: insert the approximate solution using shape functions and use the Galerkin method. Writing \(T \approx \sum_j N_j T_j\) and testing with \(N_i\) gives the weak form:

(143)\[\int_\Omega \rho c_p\, N_i N_j\, T^{n+1}_j\, d\Omega - \Delta t \int_\Omega N_i \nabla \cdot \big( k\, \nabla N_j \big)\, T^{n+1}_j\, d\Omega = \int_\Omega \rho c_p\, N_i N_j\, T^{n}_j\, d\Omega\ \ \ \ \ \ \ i=1,2,...,n\]

we proceed with integrating the diffusion term by parts (moving derivatives from \(T\) to the test function) and collecting the Neumann boundary term on \(\Gamma_N\).

(144)\[\int_\Omega \rho c_p\, N_i N_j\, T^{n+1}_j\, d\Omega + \Delta t \int_\Omega \nabla N_i \cdot \big( k\, \nabla N_j \big)\, T^{n+1}_j\, d\Omega = \int_\Omega \rho c_p\, N_i N_j\, T^{n}_j\, d\Omega - \Delta t \oint_{\Gamma_N} N_i\, \vec{q}\cdot\vec{n}\, d\Gamma\ \ \ \ \ \ \ i=1,2,...,n\]

We can proceed and write everything in terms of matrices:

(145)\[\begin{split}\left( \rho c_p M + \Delta t A \right ) T^{n+1} = \rho c_p M T^{n} + BC \\\end{split}\]

With the matrices defined as:

(146)\[\begin{split}\begin{align} \begin{split} M &= \int_\Omega N_i N_j\, d\Omega \ \ \ \ \ \ \ i,j=1,2,...,n\\ A &= \int_\Omega \left ( \frac{\partial N_i}{\partial x}\,k\,\frac{\partial N_j }{\partial x} + \frac{\partial N_i}{\partial y}\,k\,\frac{\partial N_j }{\partial y} \right ) d\Omega\ \ \ \ \ \ \ i,j=1,2,...,n\\ \end{split} \end{align}\end{split}\]

The matrix \(M\) is called the mass matrix. The terms in brackets on the LHS of equation (145) will become the new matrix that is assembled per element and added to the global stiffness matrix. If \(\rho c_p\) or \(k\) vary spatially, keep them inside the integrals (pulling them out assumes element-wise constants).

Implementation

We implement the transient behavior into our triangle script from the previous lecture. If you didn’t complete it, you can download it from here (2d_fem_steady_state_triangle.py) and the shape function definitions from here (shapes_tri.py).

We will have to make several changes to the code:

Time loop and output writing

For transient problems, we will need a time loop over all time steps and a way to store/visualize the evolving solution. So far, we have plotted only steady-state solutions that we could directly plot using matplotlib. Now we will have to come up with a different strategy as we want to visualize and analyze the complete transient solution. One good way is to use the meshio python package to save the solution in XDMF/HDF5 (VTK) format and then analyze it using ParaView. If you don’t have ParaView installed, now is the time ;)

Load the transient.xmf file in ParaView (with the Xdmf3 Reader S) and use the play button to see the transient evolution of the temperature field.

Here is a minimal code snippet that illustrates the logic:

Listing 2 Output writing.

import meshio

# mesh information (x,y,0)
points = np.hstack((GCOORD, np.zeros((GCOORD.shape[0], 1))))
cells = [("triangle", EL2NOD)]

dt = 0.025
nt = 40

with meshio.xdmf.TimeSeriesWriter('transient.xmf') as writer:
    writer.write_points_cells(points, cells)

    # model time loop
    for it in range(nt):
        # ... FEM update here ...

        # example cell data
        U = np.c_[Q_x, Q_y, np.zeros_like(Q_x)]
        writer.write_data(it*dt, point_data={"T": T}, cell_data={"U": [U], "K": [Kel]})

Note the time loop over all time steps. The number of time steps and the time step itself are chosen to make the final results look nice—they don’t have a physical meaning for the time being.

Note also, that we don’t need the python plotting at the end of the script anymore.

Matrix assembly

If we look at equation (146), we notice that we have to change the matrix assembly to 1) account for the mass matrix in the element stiffness matrix, and 2) to integrate the old temperatures into the force vector. This can be done like this:

Listing 3 matrix assembly.

# 4. compute element matrix (mass + diffusion)
# mass matrix
Me = np.outer(N, N)
# diffusion stiffness matrix
Ke = dNdx.T @ dNdx
# assemble element matrix
Ael += (rho*cp*Me + dt*Kel[iel]*Ke) * detJ * weights[ip]

# 5. assemble right-hand side from previous time step
Rhs_el += rho*cp * (Me @ T[EL2NOD[iel, :]]) * detJ * weights[ip]

Notice how the logic for the element thermal conductivity has changed - and that we need two additional physical parameters \(\rho\) and \(c_p\) .

Listing 4 model parameters.

rho         = 1
cp          = 1

Kel    = np.ones(nel)*k1
Kel[np.where(Phases==100)] = k2

Make sure that the new logic for Kel is also used in the post-processing step when computing heat fluxes.

Remark (mass matrix): The expressions above use the consistent mass \(M_{ij}=\int N_i N_j\,d\Omega\). A lumped (diagonal) mass is also common and can simplify time stepping.

Remark (boundary conditions): Neumann heat fluxes contribute to the right-hand side via the boundary integral; Dirichlet temperatures are imposed as usual by modifying rows/columns.

Results of transient diffusion problem.