# Two-field mixed $hp$-finite elements for time-dependent problems
in the refined theories of thermodynamics

Balázs Tóth  Institute of Applied Mechanics, University of Miskolc, Miskolc-
Egyetemváros, H-3515, Hungary<EMAIL_ADDRESS>Zsombor Molnár  Institute
of Applied Mechanics, University of Miskolc, Miskolc-Egyetemváros, H-3515,
Hungary<EMAIL_ADDRESS>Róbert Kovács  Department of
Energy Engineering, Faculty of Mechanical Engineering, Budapest University of
Technology and Economics, Műegyetem rkp. 3., H-1111, Hungary
Department of Theoretical Physics, Wigner Research Centre for Physics,
Budapest, Hungary<EMAIL_ADDRESS>
###### Abstract

Thanks to modern manufacturing technologies, heterogeneous materials with
complex inner structures (e.g., foams) can be easily produced. However, their
utilization is not straightforward, as the classical constitutive laws are not
necessarily valid. According to various experimental observations, the
Guyer–Krumhansl equation stands as a promising candidate to model such complex
structures. However, the practical applications need a reliable and efficient
algorithm that is capable of handling both complex geometries and advanced
heat equations. In the present paper, we present the development of a
$hp$-type finite element technique, which can be reliably applied. We
investigate its convergence properties for various situations, being
challenging in relation to stability and the treatment of fast propagation
speeds. That algorithm is also proved to be outstandingly efficient, providing
solutions four magnitudes faster than commercial algorithms.

###### keywords:

Two-field variational formulations, $hp$-version FEMs,
Maxwell–Cattaneo–Vernotte heat equation, Guyer–Krumhansl heat equation,
Transient analyzes

††journal: International Journal of Engineering Science

##  1 Introduction

The engineering practice utilizes several models to describe the material
behaviour, these equations are called constitutive equations. The most
frequent are the Navier–Stokes, Fourier and Hooke laws for continuum objects.
The common point among them is that all define an equality:

$\displaystyle
q=-\lambda\partial_{x}T,\quad\Pi=-\nu\partial_{x}v,\quad\sigma=E\varepsilon,$
(1)

in which $q,\lambda$ and $T$ are the heat flux, thermal conductivity and
temperature field; $\Pi$ describes the dynamic pressure as a consequences of
the presence of velocity gradient $\partial_{x}v$; and the stress tensor
$\sigma$ is proportional with the deformation $\varepsilon$. This simple
structure usually allows to easily eliminate the current densities and
prescribing the boundary conditions (BCs) using only one field variable, for
instance, the temperature for heat conduction problems. It has also eased the
development and implementation of various numerical and analytical solution
techniques.

However, for advanced transport models, the situation becomes more difficult.
The constitutive equation is not an equality but becomes a partial
differential equation. Therefore, the conventional initial and BCs do not work
without any further conditions [1, 2], and a more careful treatment is
necessary. A glaring example is related to the so-called Guyer–Krumhansl (GK)
equation,

$\tau\,\dot{q}+q+\lambda\,T^{\prime}-\kappa^{2}\,q^{\prime\prime}=0\quad\mathrm{in}\;\Omega=[x|x\in(0,\ell)]$
(2)

in which the partial time derivative of the heat flux ($\dot{q}$) appears with
a coefficient $\tau$ called relaxation time. Moreover, Eq. (2) consists of the
second-order spatial derivative of $q$ as well, denoted by $q^{\prime\prime}$.
Although the GK equation is usually known from the kinetic theory describing
wave-type heat conduction phenomena in crystals [3], it also can be derived on
a continuum basis [4]. Thus the coefficient $\kappa^{2}$ is not related to the
mean free path of phonons as in that continuum approach, the resulting model
is independent of the micro-scale heat transfer mechanisms [4]. Adding that
the GK equation is successfully tested in heat conduction experiments on
various rocks and foams as challenging problems [5, 6, 7], this allows us to
consider the GK equation as a promising candidate beyond Fourier and use it
for more practical engineering applications. Hence it is essential to
understand the properties of the GK equation and develop efficient and
reliable numerical solution techniques.

On that basis, a finite difference approach has been elaborated recently [8],
but it is strongly limited to simple (regular) geometries. That numerical
methodology is also adapted for rheological and nonlinear models [9], with
carefully investigating the role of initial and boundary conditions
analytically, too [1, 2]. Therefore, we aim to develop a $hp$-version finite
element method (FEM), which is able to preserve the advantageous properties of
the earlier finite difference scheme to handle the initial and boundary
conditions reliably, and, also, to be adaptable for complex geometries and
possess high accuracy and fast convergence properties.

The simplest classical FE techniques use linear piecewise polynomials to
approximate the solution and achieve the desired accuracy with mesh
refinement. The philosophy of using low-order polynomials over successively
finer meshes is called $h$-type approximation technique [10, 11, 12, 13]. In
most cases [14, 15, 16, 17], that $h$-type approach is utilised for the
Maxwell–Cattaneo–Vernotte- (MCV, with $\kappa^{2}=0$ in Eq. (2)), or for the
dual-phase-lag (DPL) equations, which models have much less practical
interest.

However, the conventional $h$-FEMs can provide slow convergences and low
accuracy for specific types of initial and BCs, as well as when certain
material and/or geometric parameters of the considered model problem is close
to their limit value. In order to overcome these numerical difficulties, the
$p$-version FE technology will be used as alternative, promising strategy,
being originally-introduced in [18, 19]. The idea is to keep the coarse mesh
fixed, and the convergence and the high accuracy are achieved solely by
increasing the polynomial degree $p$ of the approximated variables. It was
proven that the rate of convergence is much higher for $p$-FEM than that is
possible with $h$-FEM, and exponential for smooth solutions and even for non-
smooth solutions using properly chosen mesh refinement coupled with the
increase of the polynomial degree $p$ ($hp$-strategy) [19]. This is the main
motivation behind our research, thus the $p$-type approximation technique is
chosen to develop new FEM for the above-mentioned, refined theories of
thermodynamics.

### Problem formulation

The GK equation (2) is accompanied by the balance of internal energy for heat
conduction processes, that is,

$\rho\,c_{V}\,\dot{T}+q^{\prime}=0\quad\mathrm{in}\;\Omega=[x|x\in(0,\ell)],$
(3)

where $\rho$ and $c_{V}$ are the mass density and specific heat; the source
terms are neglected. In the following, we present the discretization procedure
for that system of equations.

For the sake of generality the system of basic differential Eqs. (2)–(3) is
subjected the spatial descriptions

$\displaystyle T(t,0)=\tilde{T}(t)\>,$ (4) $\displaystyle
q(t,\ell)=\tilde{q}(t)$ (5)

as Dirichlet- and Neumann-type BCs at $x=0$ and $x=\ell$, respectively, as
well as the temporal prescriptions

$T(0,x)=T_{0}(x)\quad\mathrm{and}\quad q(0,x)=q_{0}(x)$ (6)

as initial conditions at $t_{0}=0$ s.

##  2 Two-field mixed $hp$-version finite element method

In this section, new $hp$-FEMs will be presented for the numerical solution of
the MCV- and the GK model which are based on two-field mixed variational
formulations, see some similar procedures in [20, 21, 22].

###  2.1 Variational formulations

####  2.1.1 Maxwell–Cattaneo–Vernotte model

In this model problem, after testing Eq. (2), i.e., multiplying this with the
test function $v$ and integrating it over the domain $\Omega$ with setting
$\kappa^{2}=0$, then relaxing Eq. (3), i.e., multiplying this with the test
function $u$ and integrating it by parts, as well as building-in the BC (5)
and the homogeneous form of the BC (4) for $u$, we arrived at a two-field
mixed variational formulation which reads as follows. Find the duet $T\in
H^{1}(\Omega)$ and $q\in L^{2}(\Omega)$ as trial functions satisfying a priori
the BC (4) and the initial condition (6) such that

$\displaystyle-\int_{\Omega}\rho\,c_{V}\,\dot{T}\,u\,\mathrm{d}\Omega+\int_{\Omega}q\,u^{\prime}\,\mathrm{d}\Omega$
$\displaystyle=\tilde{q}_{\ell}\,u(\ell)\quad\forall u\in H^{1}(\Omega)\>,$
(7)
$\displaystyle\int_{\Omega}\tau\,\dot{q}\,v\,\mathrm{d}\Omega+\int_{\Omega}q\,v\,\mathrm{d}\Omega+\int_{\Omega}\lambda\,T^{\prime}\,v\,\mathrm{d}\Omega$
$\displaystyle=0\quad\forall v\in L^{2}(\Omega)$ (8)

satisfying a priori the homogeneous form of the BC (4). Here, $H^{1}(\Omega)$
stands for the Sobolev space of order 1 [23] and represents the regularity
property for $T$ and $u$ while $L^{2}(\Omega)$ defines square integrable
function space for $q$ and $v$.

####  2.1.2 Guyer–Krumhansl model

In this model problem relaxing again Eq. (3) with the test function $u$, but
now, after having tested Eq. (2) with test function $v$, integrating-by-parts
the last term of this variational equation, we obtain again a two-field mixed
variational formulation but with a bit different regularity property from Eq.
(8). Accordingly, we seek the variable pair $(T,\,q)\in H^{1}(\Omega)$
satisfying a priori the BC (4) and the initial condition (6) in such a way
that

$\displaystyle-\int_{\Omega}\rho\,c_{V}\,\dot{T}\,u\,\mathrm{d}\Omega+\int_{\Omega}q\,u^{\prime}\,\mathrm{d}\Omega$
$\displaystyle=\tilde{q}_{\ell}\,u(\ell)\quad\forall u\in H^{1}(\Omega)\>,$
(9)
$\displaystyle\int_{\Omega}\tau\,\dot{q}\,v\,\mathrm{d}\Omega+\int_{\Omega}q\,v\,\mathrm{d}\Omega+\int_{\Omega}\lambda\,T^{\prime}\,v\mathrm{d}\Omega+\int_{\Omega}\kappa^{2}\,q^{\prime}\,v^{\prime}\,\mathrm{d}\Omega$
$\displaystyle=q^{\prime}(\ell)\,v(\ell)-q^{\prime}(0)\,v(0)\quad\forall v\in
H^{1}(\Omega)$ (10)

hold true, ensuring a priori the homogeneous form of the BC (4) once again.
Notably, again, for both variables $(T,\,q)\in H^{1}(\Omega)$, which stands as
a crucial difference for other FE approaches for heat equations.

###  2.2 $hp$-finite element discretizations

Let us consider now the $hp$-FE discretization of the domain $\Omega$.
Thenceforward $\Omega$ is divided into $n$ physical sub-domain. Then the
master element $\Omega_{\mathrm{mas}}:=\\{\,\eta\,|\,\eta\in(-1,1)\\}$ is
mapped onto the $i$-th physical element
$K:=\\{x\,|x^{i}\in(x_{i},x_{i+1})\\}\subset\Omega_{\mathrm{h}}$ with the
nodal points $x_{i}$ and $x_{i+1}$ of the FE mesh
$\Omega_{\mathrm{h}}:0<x_{1}<x_{2}<...<x_{i}<x_{i+1}<...<x_{n}<x_{n+1}=\ell$
by the mapping function $x^{i}:=N_{1}(\eta)x_{i}+N_{2}(\eta)x_{i+1}$, where
$N_{1}=(1-\eta)/2$ and $N_{2}=(1+\eta)/2$, $i=1,\ldots,n$.

The approximation spaces for the trial- and test function group, ($T,\,q$) and
($u,\,v$) are $hp$-FE function spaces consisting of piece-wise continuous
polynomial functions, being spanned over the $i$-th element by the external
shape functions $N_{1}$ and $N_{2}$ for $p=1$, as well as their supplemented
set with the bubble modes
$N_{k}(\eta)=[L_{k-1}(\eta)-L_{k-3}(\eta)]/\sqrt{2(2k-3)}$ for $p\geq 2$,
where $L_{k}(\eta)$ are the orthogonal Legendre polynomials,
$k=3,4,\ldots,p+1$ [18, 24, 25]. The independent trial- and test functions,
$T,\;u$ and $q,\;v$ are approximated on the $i$-th element of
$\Omega_{\mathrm{h}}$ by the polynomial degrees $p+1$ and $p$, or $p+1$,
respectively, for the MCV- and the GK model according to Table 1, keeping in
mind the regularity assumptions on the trial functions $T$ and $q$, as well as
the test functions $u$ and $v$ in Eqs. (7)–(8) for the MCV model and in Eqs.
(9)–(10) for the GK model. Here $p$ is the actual polynomial degree set to
each element.

Model Function | $T,\,u$ | $q,\,v$
---|---|---
MCV | $p+1$ | $p$
GK | $p+1$ | $p+1$
Table 1: Polynomial approximation spaces.

Let us represent now the variational equations (7)–(8) and (9)–(10) in matrix
form. For the sake of simplicity let ${}^{\mathrm{h}}\boldsymbol{t}$ and
${}^{\mathrm{h}}\boldsymbol{q}$ indicate the column matrices corresponding to
the time-dependent unknown coefficients of the temperature and the heat flow,
respectively. After carrying out the numerical integrations on each element
$i$, on the spatial assembling process we have the following time-dependent
matrix equation system

$\begin{array}[]{rl}\mathcal{C}\>^{\mathrm{h}}\dot{\boldsymbol{t}}+\mathcal{Q}^{T}\>{}^{\mathrm{h}}\boldsymbol{q}&=\boldsymbol{d}\>,\\\
\mathcal{T}\>^{\mathrm{h}}\dot{\boldsymbol{q}}+\mathcal{Q}\>^{\mathrm{h}}\boldsymbol{t}+\mathcal{K}\>^{\mathrm{h}}\boldsymbol{q}&=\boldsymbol{0}\>,\end{array}$
(11)

which can be written in the simplified form

$\mathcal{A}\>^{\mathrm{h}}\dot{\boldsymbol{\alpha}}+\mathcal{B}\>^{\mathrm{h}}\boldsymbol{\alpha}=\boldsymbol{f}\>,$
(12)

where

${}^{\mathrm{h}}\boldsymbol{\alpha}=\left[\begin{array}[]{c}{}^{\mathrm{h}}\boldsymbol{t}\\\
{}^{\mathrm{h}}\boldsymbol{q}\end{array}\right]\>,\quad\mathcal{A}=\left[\begin{array}[]{cc}\mathcal{C}&\boldsymbol{0}\\\
\boldsymbol{0}&\mathcal{T}\end{array}\right]\;,\quad\mathcal{B}=\left[\begin{array}[]{cc}\boldsymbol{0}&\mathcal{Q}^{T}\\\
\mathcal{Q}&\mathcal{K}\end{array}\right]\quad\mathrm{and}\quad\boldsymbol{f}=\left[\begin{array}[]{c}\boldsymbol{d}\\\
\boldsymbol{0}\end{array}\right]\;,$ (13)

in which $\mathcal{C}$, $\mathcal{T}$ and $\mathcal{K}$ denote the consistent
matrices of the specific heat, the relaxation term and the heat conductivity,
respectively, while $\mathcal{Q}$ is the consistent coupling matrix of the
system.

##  3 Numerical experiments

In this section the newly-developed two-field $hp$-version FEMs will be tested
on a representative initial-boundary value problem, focusing on the transient
analyzes. In order to demonstrate the computational performance of the
$hp$-type mixed FEs, first we investigate the point-wise $h$\- and
$p$-convergence of the relative error that are measured, respectively, in the
maximum norm of the temperature and the heat flow

$e^{\mathrm{hp}}=\frac{\underset{t\in[t_{0},t_{1}]}{\max}\left|\>{}^{\mathrm{h}}\boldsymbol{t}-\boldsymbol{t}_{\mathrm{ref}}\right|}{\underset{t\in[t_{0},t_{1}]}{\max}\left|\boldsymbol{t}_{\mathrm{ref}}\right|}\quad\textrm{and}\quad
e^{\mathrm{hp}}=\frac{\underset{t\in[t_{0},t_{1}]}{\max}\left|\>{}^{\mathrm{h}}\boldsymbol{q}-\boldsymbol{q}_{\mathrm{ref}}\right|}{\underset{t\in[t_{0},t_{1}]}{\max}\left|\boldsymbol{q}_{\mathrm{ref}}\right|}\;.$
(14)

The related thermodynamic system (12) is solved numerically for the relatively
small time interval $t\in[0,10]$ s, using the implicit time integration scheme
[12]. The number of the constant time step is set to $n_{t}=10000$, yielding
the time step size $\Delta t=0.001$ s. The considered bodies of length
$\ell=0.005$ m are made of rock-like materials which have
$c_{V}=800\;\mathrm{(J/kgK)}$, $\rho=2600\;\mathrm{kg/m^{3}}$. Additionally,
the new parameters $\tau$ and $\kappa^{2}$ are determined based on earlier
experimental data [5] to be realistic, therefore we choose
$\tau=\\{0.05,0.15,0.3\\}$ s, with $\kappa^{2}=8\cdot 10^{-6}$ m2.

However, in order to challenge the numerical procedure, we also test for
unrealistically high $\kappa^{2}=0.8$ m2 as that coefficient greatly
influences the characteristic speed of diffusion and makes the propagation
faster with five magnitudes with respect to the previous situations.

Furthermore, these are subjected to the time-dependent Neumann-type BCs

$\tilde{q}_{0}(t)=10000\,\frac{c_{1}\,c_{2}}{c_{2}-c_{1}}\left[\exp\left(-c_{1}\frac{t}{t_{p}}\right)-\exp\left(-c_{2}\frac{t}{t_{p}}\right)\right]\quad\mathrm{and}\quad\tilde{q}_{\ell}(t)=0$
(15)

as prescribed heat flows and the initial conditions $T_{0}(x)=293$ K and
$q_{0}(x)=0$, where $c_{1}=1/0.075$, $c_{2}=6$, $t_{p}=0.008$ s, namely, these
are initially at rest, i.e, in equilibrium state. The $c_{1}$ and $c_{2}$
coefficients are chosen based on [6], to keep the excitation to be realistic
as much as possible.

When the computational performance of the $p$-approximation is analyzed, a 8-,
52- and 20-element uniform mesh kept fixed on the domain $\Omega$ for the GK
model with the parameter settings $\kappa^{2}=0.8$ m2, $\kappa^{2}=0.000008$
m2 and the MCV model, respectively, while the polynomial degree $p$ is ranging
from 2 to 8 for each element. During the $h$-convergence studies, the domain
$\Omega$ is uniformly-refined in 7 steps from the element number $n=8$, 52 and
20 to 20, 88 and 44, respectively, for the GK model with $\kappa^{2}=0.8$ m2,
$\kappa^{2}=0.000008$ m2 and the MCV model. In each step the polynomial degree
$p$ is set to 2 and kept fixed on each element. At each step, the number of
degrees of freedom (DOF) is calculated as the sum of the total unknown
coefficients occurring in $\boldsymbol{\alpha}$.

The relative error-convergence curves of the front- and rear side temperature,
$T(t,0)$ and $T(t,\ell)$, as well as the mid-side heat flow $q(t,\ell/2)$
obtained for the MCV model and the GK model with relatively large and small
value of $\kappa^{2}$ (0.8 m2 and 0.000008 m2) are plotted against the number
of DOF on log-log scales in Figs. 1 and 2–3 for the relaxation time $\tau=0.3$
s, $\tau=0.15$ s and $\tau=0.05$ s, respectively.

Following from the characteristics of the relative error-curves, monotonically
high rates are experienced both for $h$\- and for $p$-convergence. As
expected, the $p$-convergence is much faster than the $h$-convergence. Namely,
the higher-order polynomial approximation helps a lot to increase the
convergence rate. Thus, the desired accuracy is achieved with the use of much
less number of DOF by the $p$-approximation. In the asymptotic range, the
exponential type of the $p$-convergence behavior is observed while the
$h$-convergence shows rather an algebraic type of convergence.

It can also be experienced that there is no significant influence of the
relaxation time $\tau$ on the convergence rates. However, the convergence
curves are shifted down (without changing their slopes) as $\tau$ is decreased
achieving a lower relative error value. Besides, it can be seen that the
computational limit is reached at a very small error level (at approximately
$10^{-8}$).

Giving now as the illustration of the $hp$-FE solutions for the MCV model and
the GK model with the quite small value of $\kappa^{2}$ (0.000008 m2), the
histories of the dimensionless rear- and front side temperatures are depicted,
separately for three different relaxation times $\tau=0.3$ s, $\tau=0.15$ s
and $\tau=0.05$ s. There is no any computational reason behind the
dimensionless rescaling of temperature, this is used only for better
overviewing the plots as the model is linear, and solely the characteristics
are important. The adiabatic steady-state temperature is set to $1$ for these
Figures. In Figs. 4–6 and, 8–10, respectively, using a 100-element-mesh with
the relatively high polynomial degree $p=10$ and comparing all the results to
the $hp$-FE solution of the classical, Fourier model. Figs. 7 and 11 represent
the time series of the dimensionless rear- and front side temperature for
over-diffuse thermodynamic system, i.e., for $\kappa^{2}=0.8$ m2, zooming the
$hp$-FE solutions in the rapidly change and very small time region $[0,0.02]$
s. These figures exhibit very well that the $\tau$-value has no significant
effect on the solution for over-diffuse system.

$10^{2}$$10^{2.5}$$10^{3}$$10^{-3}$$10^{-2}$number of DOFrelative error in maximum normFront side temperature | $10^{2}$$10^{2.5}$$10^{3}$$10^{-11}$$10^{-9}$$10^{-7}$$10^{-5}$number of DOFRear side temperature | $10^{2}$$10^{2.5}$$10^{3}$$10^{-5}$$10^{-3}$$10^{-1}$number of DOFMid-side heat flow$\tau=0.3\;\mathrm{s},\;h$-conv.$\tau=0.15\;\mathrm{s},\;h$-conv.$\tau=0.05\;\mathrm{s},\;h$-conv.$\tau=0.3\;\mathrm{s},\;p$-conv.$\tau=0.15\;\mathrm{s},\;p$-conv.$\tau=0.05\;\mathrm{s},\;p$-conv.
---|---|---
Figure 1: Convergence histories of the relative errors measured in maximum norm for the front- and rear side temperature, as well as the mid-side heat flow – MCV model $10^{1.5}$$10^{2}$$10^{-9}$$10^{-8}$$10^{-7}$$10^{-6}$$10^{-5}$number of DOFrelative error in maximum normFront side temperature | $10^{1.5}$$10^{2}$$10^{-10}$$10^{-9}$$10^{-8}$$10^{-7}$$10^{-6}$number of DOFRear side temperature | $10^{1.5}$$10^{2}$$10^{-8}$$10^{-6}$$10^{-4}$$10^{-2}$number of DOFMid-side heat flow$\tau=0.3\;\mathrm{s},\;h$-conv.$\tau=0.15\;\mathrm{s},\;h$-conv.$\tau=0.05\;\mathrm{s},\;h$-conv.$\tau=0.3\;\mathrm{s},\;p$-conv.$\tau=0.15\;\mathrm{s},\;p$-conv.$\tau=0.05\;\mathrm{s},\;p$-conv.
---|---|---
Figure 2: Convergence histories of the relative errors measured in maximum norm for the front- and rear side temperature, as well as the mid-side heat flow – GK model with $\kappa^{2}=0.8\;\mathrm{m}^{2}$ $10^{2.5}$$10^{3}$$10^{-10}$$10^{-8}$$10^{-6}$$10^{-4}$$10^{-2}$number of DOFrelative error in maximum normFront side temperature | $10^{2.5}$$10^{3}$$10^{-10}$$10^{-9}$number of DOFRear side temperature | $10^{2.5}$$10^{3}$$10^{-9}$$10^{-7}$$10^{-5}$$10^{-3}$number of DOFMid-side heat flow$\tau=0.3\;\mathrm{s},\;h$-conv.$\tau=0.15\;\mathrm{s},\;h$-conv.$\tau=0.05\;\mathrm{s},\;h$-conv.$\tau=0.3\;\mathrm{s},\;p$-conv.$\tau=0.15\;\mathrm{s},\;p$-conv.$\tau=0.05\;\mathrm{s},\;p$-conv.
---|---|---
Figure 3: Convergence histories of the relative errors measured in maximum
norm for the front- and rear side temperature, as well as the mid-side heat
flow – GK model with $\kappa^{2}=0.000008\;\mathrm{m}^{2}$ Figure 4: Time
series of the dimensionless rear side temperature for the MCV- and GK model
with $\kappa^{2}=0.000008\;\mathrm{m}^{2}$ – $\tau=0.3$ s, $p=10$ and $n=100$
Figure 5: Time series of the dimensionless rear side temperature for the MCV-
and GK model with $\kappa^{2}=0.000008\;\mathrm{m}^{2}$ – $\tau=0.15$ s,
$p=10$ and $n=100$ Figure 6: Time series of the dimensionless rear side
temperature for the MCV- and GK model with
$\kappa^{2}=0.000008\;\mathrm{m}^{2}$ – $\tau=0.05$ s, $p=10$ and $n=100$
Figure 7: Influence of $\tau$ on the time series of the dimensionless rear
side temperature for the GK model with the over-diffuse setting
$\kappa^{2}=0.8\;\mathrm{m}^{2}$ – $p=10$ and $n=100$ Figure 8: Time series of
the dimensionless front side temperature for the MCV- and GK model with
$\kappa^{2}=0.000008\;\mathrm{m}^{2}$ – $\tau=0.3$ s, $p=10$ and $n=100$
Figure 9: Time series of the dimensionless front side temperature for the MCV-
and GK model with $\kappa^{2}=0.000008\;\mathrm{m}^{2}$ – $\tau=0.15$ s,
$p=10$ and $n=100$ Figure 10: Time series of the dimensionless front side
temperature for the MCV- and GK model with
$\kappa^{2}=0.000008\;\mathrm{m}^{2}$ – $\tau=0.05$ s, $p=10$ and $n=100$
Figure 11: Influence of $\tau$ on the time series of the dimensionless front
side temperature for the GK model with the over-diffuse setting
$\kappa^{2}=0.8\;\mathrm{m}^{2}$ – $p=10$ and $n=100$

##  4 Discussion

Here we presented a novel numerical approach, specifically developed for the
Guyer–Krumhansl equation as that model could have significant practical
interest in engineering based on the available experimental data. We
successfully demonstrated that the solutions converge, thus are stable and
consistent and reproduce the experimentally observed required characteristics,
on contrary to COMSOL [8]. That is essential, since the usual approaches do
not work, it is not possible to use the temperature as a single field variable
for heat flux BCs [2].

Moreover, that technique is incredibly more efficient than the one offered by
COMSOL, runs about 4 magnitudes faster thanks to the low number of DOFs and
elements. That advantage becomes more significant in two and three-dimensional
problems, therefore our next aim is to improve and implement the present
approach for complex geometries in higher spatial dimensions. Additionally,
due to the global energy crisis and chip shortage, the efficiency of
algorithms has increasing importance. That huge increase in speed would enable
the real-time monitoring of heterogeneous materials with complex inner
structure without massive computational capacity.

##  5 Acknowledgements

The research reported in this paper is part of project no. BME-NVA-02,
implemented with the support provided by the Ministry of Innovation and
Technology of Hungary from the National Research, Development and Innovation
Fund, financed under the TKP2021 funding scheme. This paper was supported by
the János Bolyai Research Scholarship of the Hungarian Academy of Sciences.
The research reported in this paper and carried out at BME has been supported
by the grants National Research, Development and Innovation Office-NKFIH FK
134277.

## References

  * Kovács [2018] R. Kovács,  Analytic solution of Guyer-Krumhansl equation for laser flash experiments,  International Journal of Heat and Mass Transfer 127 (2018) 631–636.
  * Kovács [2022] R. Kovács,  Analytical treatment of nonhomogeneous initial states for non-Fourier heat equations,  International Communications in Heat and Mass Transfer 134 (2022) 106021\.
  * Guyer and Krumhansl [1966] R. A. Guyer, J. A. Krumhansl,  Thermal Conductivity, Second Sound, and Phonon Hydrodynamic Phenomena in Nonmetallic Crystals,  Physical Review 148 (1966) 778–788.
  * Ván and Fülöp [2012] P. Ván, T. Fülöp,  Universality in heat conduction theory – weakly nonlocal thermodynamics,  Annalen der Physik (Berlin) 524 (2012) 470–478.
  * Both et al. [2016] S. Both, B. Czél, T. Fülöp, G. Gróf, Á. Gyenis, R. Kovács, P. Ván, J. Verhás,  Deviation from the Fourier law in room-temperature heat pulse experiments,  Journal of Non-Equilibrium Thermodynamics 41 (2016) 41–48.
  * Fehér and Kovács [2021] A. Fehér, R. Kovács,  On the evaluation of non-Fourier effects in heat pulse experiments,  International Journal of Engineering Science 169 (2021) 103577. ArXiv:2101.01123.
  * A.Lunev et al. [2022] A.Lunev, A. Lauerer, V. Zborovskii, F. Léonard,  Digital twin of a laser flash experiment helps to assess the thermal performance of metal foams,  International Journal of Thermal Sciences 181 (2022) 107743.
  * Rieth et al. [2018] A. Rieth, R. Kovács, T. Fülöp,  Implicit numerical schemes for generalized heat conduction equations,  International Journal of Heat and Mass Transfer 126 (2018) 1177 – 1182.
  * Pozsar et al. [2020] A. Pozsar, M. Szücs, R. Kovács, T. Fülöp,  Four spacetime dimensional simulation of rheological waves in solids and the merits of thermodynamics,  Entropy 22 (2020) 1376.
  * Bathe [1996] K.-J. Bathe, Finite Element Procedures, Prentice Hall, Upper Saddle River, New Jersey, 1996.
  * Ciarlet [1978] P. G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland Publishing Company, Amsterdam, 1978.
  * Hughes [1987] T. J. R. Hughes, The Finite Element Method, Linear Static and Dynamic Finite Element Analysis, Prentice Hall, Englewood Cliffs, New Jersey, 1987.
  * Reddy [2002] J. N. Reddy, Energy Principles and Variational Methods in Applied Mechanics, 2nd ed., John Wiley & Sons, New York, 2002\.
  * Deka and Dutta [2020] B. Deka, J. Dutta,  Finite element methods for non-fourier thermal wave model of bio heat transfer with an interface,  Journal of Applied Mathematics and Computing 62 (2020) 701–724.
  * Yang et al. [2019] Z.-B. Yang, Z.-K. Wang, S.-H. Tian, X.-F. Chen,  Analysis and modelling of non-fourier heat behavior using the wavelet finite element method,  Materials 12 (2019).
  * Nazmdeh et al. [2021] H. Nazmdeh, M. Vahabi, M. A. Nazari,  Finite element modeling of Non-Fourier heat transfer in a cancerous tissue with an injected fat layer during hyperthermia treatment,  Journal of Thermal Biology 100 (2021) 103073.
  * Xu and Li [2003] B. Xu, B. Q. Li,  Finite element solution of non-Fourier thermal wave problems,  Numerical Heat Transfer, Part B: Fundamentals 44 (2003) 45–60.
  * Szabó and Babuška [2011] B. Szabó, I. Babuška, Introduction to Finite Element Analysis: Formulation, Verification and Validation, John Wiley & Sons, New York, 2011.
  * Szabó et al. [1981] B. Szabó, I. Babuška, I. N. Katz,  The $p$-version of the finite element method,  SIAM Journal on Numerical Analysis 18 (1981) 515–545.
  * Tóth [2018] B. Tóth,  Dual and mixed nonsymmetric stress-based variational formulations for coupled thermoelastodynamics with second sound effect,  Continuum Mechanics and Thermodynamics 30 (2018) 319–345.
  * Tóth [2016] B. Tóth,  Multi-field dual-mixed variational principles using non-symmetric stress field in linear elastodynamics,  Journal of Elasticity 122 (2016) 113–130.
  * Tóth [2021] B. Tóth,  Natural frequency analysis of shells of revolution based on hybrid dual-mixed $hp$-finite element formulation,  Applied Mathematical Modelling 98 (2021) 722–746.
  * Brenner and Scott [2008] S. C. Brenner, L. R. Scott, The Mathematical Theory of Finite Element Methods, 3rd ed., Springer-Verlag, New York, 2008.
  * Demkowicz [2007] L. Demkowicz, Computing with $hp$-Adaptive Finite Elements. One- and Two-Dimensional Elliptic and Maxwell Problems, volume I of Applied Mathematics and Nonlinear Science, Chapman & Hall/CRC Press, Taylor & Francis Group, 2007.
  * Düster et al. [2017] A. Düster, E. Rank, B. Szabó, The $p$-Version of the Finite Element and Finite Cell Methods, Encyclopedia of Computational Mechanics, 2nd ed., John Wiley & Sons, Ltd, 2017, pp. 1–35. doi:10.1002/9781119176817.ecm2003g.