The present document presents the derivation of a finite element code for one dimensional transient heat conduction problems using both linear and quadratic elements. The code has been implemented in a Python package called thermesh is available for download from the author’s Github page. Three validation cases are presented at the end of this document to show that the code is implemented correctly. Further, the validation cases can be used as examples and inspiration.
In the case of 1D heat conduction, the governing partial differential equation reads:
with \(T(z,t)\) the temperature at location \(z\) and time \(t\), \(ρ c\text{p}\) the volumetric heat capacity, \(k\) the thermal conductivity, and \(\dot{Q}\) an internal source or sink, e.g. due to a phase transformation. The partial differential equation is subject to boundary conditions at the two ends of the domain, which can be of the Dirichlet type:
or of the Neumann type:
with \(T_\text{left}\) and \(T_\text{right}\) an imposed temperature and \(q_\text{left}\) and \(q_\text{right}\) an imposed heat flux density. The latter could be either directly imposed \(\hat{q}\), for example as a result of laser heating, or be the result of convection:
in which \(h\) and \(T∞\) are the heat transfer coefficient and the far field temperature, or the result of thermal radiation:
where \(ε\) is the surface emissivity and \(σ\) is Stefan’s constant.
We can approximate the solution \(T(z,t)\) of the partial differential equation using the weighted residual method:
with \(w\) a weighting function.
The domain is now split up in \(N\) smaller elements, as is indicated in Figure fig:elements which shows an element \((e)\) of length \(\ell\) that is bounded by the nodes \(i\) and \(j\). The weighted residual (Equation \ref{eq:weighted_residual}) can now be rewritten as:
We will now make sure that the weighted residual vanishes for each element. Further, we get rid of the second derivative with respect to \(z\) in the second term using integration by parts:
with \(z_i\) and \(z_j\) the element end points. Realizing that the third term represents the flux (\(q = -k ∂ T / ∂ z\)), the equation can be rewritten as:
The temperature inside an element is approximated as a linear interpolation between the bounding node temperatures as:
with \(ζ\) a local coordinate, as is illustrated in Figure fig:elements, defined as:
while the two shape functions are:
The spatial derivative of the temperature with respect to \(z\) can now be calculated as:
with:
such that:
which intuitively makes sense. Further, noting that the shape functions do not depend on time, we can rewrite the time derivative as:
Following the Galerkin method, we choose our weighting function \(w\) to be our shape functions. The equation for the weighted residual for an element (Equation \ref{eq:weighted_residual_el}) can now be rewritten as:
with \(N_k\) the two shape functions as defined in Equation \ref{eq:shape_functions}.
Before evaluating the integrals, we first rewrite the expressions from the previous section into a matrix-vector form. Starting with Equation \ref{eq:T_approx}, which can be written as:
in which:
The spatial derivative of the temperature (Equation \ref{eq:dTdz}) yields:
with:
while the time-derivative of the temperature can be rewritten as:
Further, for convenience, we will write our weighting functions as:
We can now evaluate the integrals, starting with the first term:
We can rewrite this integral in terms of \(ζ\), by making use of the derivative of \(ζ\) with respect to \(z\):
such that:
with:
In the same manner, the second term yields:
with:
The first term on the right hand side yields:
with \(\dot{Q}\) the heat source term for the element between nodes \(i\) and \(j\). The second term, with the heat flux \(q\) on the boundary, is first expanded to include both a direct heat flux \(\hat{q}\) and a flux due to convection:
which yields:
The term with the direct heat flux \(\hat{q}\) is evaluated as:
with \(\hat{q}_k\) the heat flux on the \(k\)-th node. The convective term can be accounted for using a stiffness matrix for convection:
and an additional term in the force vector:
As an example for the stiffness matrix \(\mathbf{H}\), in case of a convective boundary condition at the j-th node, where \(N_i = 0\), this term would evaluate as:
which intuitively makes sense. The force vector is now combined as:
The final element equation can now be assembled by substituting Equations \ref{eq:C}, \ref{eq:K}, \ref{eq:H} and \ref{eq:f} in Equation \ref{eq:galerkin}:
With the local damping and stiffness matrices determined for each element, we can assemble the global matrices using the node locations and element connectivity in the global system.
In the case of quadratic shape functions, the temperature inside an element is approximated as:
with reference to fig:qua_elements for the node locations. The three shape functions are now defined as:
or:
The mapping between the local coordinate \(ζ\) and the global coordinate \(x\) is achieved by:
Here, for convenience, we will consider a regular element which means that:
The Jacobian is now evaluated as:
The temperature gradient with respect to \(z\) can be written as:
with:
Now we can derive the damping matrix \(\mathbf{C}\) and the stiffness matrix \(\mathbf{K}\) in the same manner as we have done for the linear shape functions:
while the term for the internal heat source \(\dot{Q}\) is evaulated as:
The final step is to integrate the equation with respect to time. For this purpose, we will discretize the temporal variable using the so-called \(Θ\)-method:
with the subscript \(n\) denoting the time step, and where \(Θ ∈ [0, 1]\). Common values of \(Θ\) are:
Equation \ref{eq:theta} can be rearranged as:
which can be solved to obtain the temperatures for each time step.
The finite elements derived here were implemented in an object-oriented Python code called thermesh, which can be found on the author’s Github page. Please note that, with the focus on readability, the code is far from optimized. Nevertheless, I believe it should still be more than fast enough for most problems. This section presents three short validation cases to show that the code is implemented correctly. Although all presented cases use linear elements, the code was also validated for the quadratic elements. Details on material properties and dimensions used for each case can be found in the respective Python files.
Consider a domain of length \(L\) with a uniform initial temperature \(T_0\). For \(t>0\) the temperature at one end is raised to a value of \(T\text{end}\), while the other end is kept at the initial temperature:
In case the initial temperature equals 0.0 \(ˆ\)C , the analytical solution\footnote{The Mathematics of Diffusion, Crank, 1975, pp 49-50.} yields:
with \(α = k/ρ c\text{p}\) the thermal diffusivity. The left graph in Figure \ref{fig:step_compare} shows the temperature distribution at different times. Code listing \ref{lst:step} illustrates how to solve this problem using thermesh. The right graph in Figure \ref{fig:step_compare} shows the finite element solution for 10 linear elements of equal length. Good comparison is obtained between the numerical and analytical solution. The code for this comparison is available in step_change.py.
Consider a semi-infinite domain with a uniform initial temperature. The domain is subjected to a constant heat flux \(\hat{q}\) at the boundary which causes the temperature to increase:
The analytical solution\footnote{Heat Transfer, Nillis \& Klein, 2008, p 362.} for the temperature increase \(Δ T\) in the domain yields:
with \(\erfc\) the complementary error function. The left graph in \ref{fig:flux_compare} shows the temperature increase at 2, 10 and 25 seconds. Code listing \ref{lst:flux} illustrates the code to solve this problem using thermesh. Of course, one has to make sure that the domain is large enough (or the time short enough) for it to be considered a semi-infinite solid. The right graph in Figure \ref{fig:flux_compare} shows the finite element solution for 10 linear elements of equal length. Good comparison is obtained between the numerical and analytical solution. The code for this comparison is available in the Python file constant_heat_flux.py.
In the last example, we again consider a semi-infinite domain with a uniform initial temperature. The domain is now subjected to a convective boundary condition with a heat transfer coefficient \(h\) and a far field temperature \(T∞\).
The analytical solution\footnote{See footnote 2.} for the temperature increase \(Δ T\) in the domain yields:
The left graph in \ref{fig:conv_compare} shows the temperature increase at 2, 10 and 25 seconds. Code listing \ref{lst:conv} illustrates the code to solve this problem using thermesh. The right graph in Figure \ref{fig:conv_compare} shows the finite element solution for 10 linear elements of equal length. As can be seen, also here good comparison is obtained between the numerical and analytical solution. The code for this comparison is available in the Python file convective_bc.py.