Transient heat diffusion
The next step is to derive and implement the unsteady (time-dependent) heat diffusion equation and to solve an example problem for the cooling of the lithosphere. The unsteady heat diffusion equation looks like this:
Two differences to the previously considered steady 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 not constant (neglected) any more. 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. To solve it we will use a “trick” from finite differences – we will simply write the time derivative in finite differences form.
Re-arrange equation (140) so that all known temperatures \(T^n\) are on the Rhs and all unknown temperatures \(T^{n+1}\) are on the Lhs.
FEM form
Now we proceed in the usual way: insert the approximate solution using shape functions and use the Galerkin method:
and integrate the diffusion term by parts:
We can proceed and write everything in terms of matrices:
With the matrices being define as:
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.
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_transient_triangle.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 vtk format and to later analyze is using Paraview If you don’t have paraview installed, no is the time ;)
Here is a code snippit that illustrates the logic:
# mesh information
points=np.hstack((GCOORD, GCOORD[:,0].reshape(-1,1)*0)) #must have 3 components (x,y,z)
cells=[("triangle",EL2NOD)]
# write initial mesh
writer=meshio.xdmf.TimeSeriesWriter('transient.xmf')
writer.__enter__() # have to add this: import hdf5 and open file ...
writer.write_points_cells(points, cells)
dt = 0.025
nt = 40
# model time loop
for t in range(0,nt):
#our FEM code here
#save results
#cell data
U=np.hstack((Q_x.reshape(-1,1),Q_y.reshape(-1,1)))
U=np.hstack((U,U[:,0].reshape(-1,1)*0))
#save data
writer.write_data(t, point_data={"T": T},cell_data={"U": [U], "K": [Kel]})
writer.__exit__() # close file
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:
# 4. compute element stiffness matrix
Ael = Ael + (rho*cp*np.outer(N,N) + dt*Kel[iel]*np.matmul(dNdx.T, dNdx))*detJ*weights[ip] # [nnodel,1]*[1,nnodel] / weights are missing, they are 1
# 5. assemble right-hand side
Rhs_el = Rhs_el + rho*cp*np.matmul(np.outer(N,N), np.take(T, EL2NOD[iel,:], axis=0 ))*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\) .
rho = 1
cp = 1
Kel = np.ones(nel)*k1
Kel[np.where(Phases==100)] = k2
Make sure that the new logical for Kel
is also used in the post-processing step when computing heat fluxes.
Results of transient diffusion problem.