1 Introduction

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

$$\begin{aligned} q=-\lambda \partial _x T, \quad \Pi = - \nu \partial _x v, \quad \sigma = E \varepsilon , \end{aligned}$$
(1)

in which \(q, \lambda \), and T are the heat flux, thermal conductivity, and temperature field; \(\Pi \) describes the dynamic pressure as a consequence 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 us to easily eliminate the current densities and prescribe the boundary conditions (BC) 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 (PDE). Therefore, the conventional initial and BCs do not work without further conditions [1, 2], and more careful treatment is necessary. A glaring example is related to the so-called Guyer–Krumhansl (GK) equation,

$$\begin{aligned} \tau \,{\dot{q}}+q+\lambda \,T'-\kappa ^2\,q''=0\quad \textrm{in}\;I\times \Omega \end{aligned}$$
(2)

in which I and \(\Omega \) are the time and space domain (defining them after Eq. (3)) as well as the partial time derivative of the heat flux (\(\dot{q}\)) appears with a coefficient \(\tau \) called the relaxation time. Moreover, Eq. (2) consists of the second-order spatial derivative of q as well, denoted by \(q''\). While it is often claimed that the GK equation can only be well-posed with additional boundary conditions (such as adding the gradient of q as well), it is proved that the same number of BCs is enough, and well-posed with q-boundaries [3]. This is also supported by the existing analytical solutions for time-dependent BCs [1, 4] and for other BCs [5, 6].

Although the GK equation is usually known from the kinetic theory describing wave-type heat conduction phenomena in crystals [7], it also can be derived on a continuum basis [8]. Thus, the coefficient \(\kappa ^2\) is not related to the mean free path of phonons as in that continuum approach, and the resulting model is independent of the micro-scale heat transfer mechanisms [8]. Adding that the GK equation is successfully tested in heat conduction experiments on various rocks and foams as challenging problems, [4, 9, 10], and 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 [11], but it is strongly limited to simple (regular) geometries. That numerical methodology is also adapted for rheological and nonlinear models [12], 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 can 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.

Furthermore, COMSOL’s ability to solve the GK equation is not satisfactory, and this is tested by utilizing COMSOL’s Mathematics Module in version v5.3a and implementing the GK equation using the local form of PDEs [11]. Three characteristic settings were benchmarked: \(\kappa ^2=0\) leading to the MCV equation, \(\kappa ^2/\tau = \alpha \) to reproduce Fourier’s solution, and \(\kappa ^2/\tau > \alpha \) to produce over-diffusive solutions, with \(\alpha = \lambda /(\rho c_V)\) being the thermal diffusivity. In the first two situations, COMSOL can provide a physically valid temperature history, although its resource requirements—CPU and RAM usage—are rather high. However, in the last one in which \(\kappa ^2/\tau > \alpha \), the resulting temperature history was found to be far from realistic [11], false solutions can easily be found for a stable and seemingly convergent method. These solutions are independent of the applied mesh and time-stepping algorithm. Figure 1 compares the physically valid temperature and the false one from COMSOL, presenting the importance and the need for a reliable, efficient FEM.

Fig. 1
figure 1

Demonstrative solutions for the GK equation for the same set of parameters [11]: \(\alpha \)=1e-4 m\(^2\)/s, \(\tau = 0.005\) s and \(\kappa ^2\)=5e-6 m\(^2\). Upper: the dimensionless temperature history is obtained with an unconditionally stable finite difference method, validated with an analytical solution. Lower: COMSOL’s solution

Full size image

We want to underline that this is also the case for simple wave propagation, the built-in wave equation module of COMSOL, using its built-in time-stepping algorithms, provides different solutions for each algorithm, and eventually, without any supplementary solution technique independent of COMSOL, one cannot find the valid, physically admissible solution [12].

The most straightforward 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 [13,14,15,16]. In most cases [17,18,19,20,21,22,23,24], that quite well-applicable h-type approaches or other approximation schemes are utilized 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. Furthermore, the time fractional version of the MCV model is also used as a refined heat conductivity model, see the theory; and finite difference schemes, as well as spectral method elaborated on its numerical solution in [25]; and [26,27,28], as well as [29]. Besides, in paper [30] and its improvement [31], an analytical solution is derived for modeling the laser short-pulse heating process of a solid body by the use of the Laplace transformation.

However, the conventional h-FEMs can provide slow convergences and low accuracy for specific types of initial and BCs, as well as when particular material and/or geometric parameters of the considered model problem are close to their limit value, see for shell problems in [32,33,34]. In order to circumvent these computational difficulties, the p-version FE technology will be used as an alternative, promising strategy, being introduced initially in [35, 36]. 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) [36]. Besides, the so-called discontinuous Petrov–Galerkin-type hp-version FEMs are also well suitable for stable and accurate transient heat conductivity equation solutions, see in [37]. This is the primary motivation behind our research. Thus, the p-type approximation technique is chosen to develop new FEM for the above-mentioned refined theories of thermodynamics.

Accordingly, the novelty of this paper is on two main accounts: (i) derivation of new two-field variational formulations for the MCV and GK model of the refined theories of thermodynamics, as well as (ii) the development of new, advanced hp-version mixed FEMs that handle the temperature and heat flux fields as independent variables.

1.1 Problem formulation

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

$$\begin{aligned} \rho \,c_V\,{\dot{T}}+q'=0\quad \textrm{in}\;I\times \Omega \,, \end{aligned}$$
(3)

where \(I\in (t_0,t_1]\) and \(\Omega =[x|x\in (0,\ell )]\) define the considered time and space domain, as well as \(\rho \) and \(c_V\) are the mass density and specific heat; the source terms are neglected. For the sake of generality, the system of basic differential Eqs. (2, 3) is subjected to the spatial descriptions

$$\begin{aligned} T(t,0)&={\tilde{T}}(t)\quad \textrm{in}\;I\,,\end{aligned}$$
(4)
$$\begin{aligned} q(t,\ell )&={\tilde{q}}(t)\quad \textrm{in}\;I\,, \end{aligned}$$
(5)

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

$$\begin{aligned} T(0,x)=T_0(x)\quad \textrm{and}\quad q(0,x)=q_0(x)\end{aligned}$$
(6)

as initial conditions at \(t_0=0\) s, which has to be compatible with the boundary conditions (4, 5) at \(x=0\) and \(x=\ell \):

$$\begin{aligned} T_0(0)={\tilde{T}}(t_0)\quad \text {and}\quad q_0(\ell )={\tilde{q}}(t_0)\,.\end{aligned}$$
(7)

The GK-system (2, 3) is a second-order partial differential equation system, more precisely, second-order with respect to the space coordinate, but first-order in time. If \(\kappa \) (or \(\kappa ^2\)) tends (or equal to) zero, the MCV-system is obtained as an entirely first-order partial differential equation system.

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 is based on two-field mixed variational formulations

2.1 Variational formulations

In this subsection, two new variational formulations will be derived for the MCV and the GK model problem of the refined theory of thermodynamics, serving as a basis for the hp-FEMs. Some similar mathematical procedures are elaborated on mixed variational theorems of coupled time-dependent thermoelasticity problems in [32, 38].

2.1.1 Maxwell–Cattaneo–Vernotte model

In this model problem, as the first step of the derivation, after having divided the constitutive Eq. (2) by \(\lambda \) and testing it, i.e., multiplying it with the time-independent test function v(x) and integrated it over the domain \(\Omega \) with setting \(\kappa ^2=0\) we obtain its weakened form

$$\begin{aligned} \int _{\Omega }\left( \frac{\tau }{\lambda }\,{\dot{q}}+\frac{q}{\lambda }+\,T'\,\right) v\,\textrm{d}\Omega =0\,. \end{aligned}$$
(8)

In the next step of the derivation, multiplying the energy balance equation (3) with another time-independent test function u(x) and integrating it over the domain \(\Omega \) we have

$$\begin{aligned} \int _{\Omega }\left( \rho \,c_V\,{\dot{T}}+q'\,\right) u\,\textrm{d}\Omega =0\,. \end{aligned}$$
(9)

Now spatially integrating-by-parts this equation, as well as building in the BC (5) and the homogeneous form of the BC (4) for the test function u, the second term of the latter integral becomes

$$\begin{aligned} \int _{\Omega }q'\,u\,\textrm{d}\Omega =\left[ q\,u\right] _0^{\ell }-\int _{\Omega }q\,u'\,\textrm{d}\Omega ={\tilde{q}}_{\ell }\,u(\ell )-\int _{\Omega }q\,u'\,\textrm{d}\Omega \,. \end{aligned}$$
(10)

This process is called relaxation. Accordingly, upon substitution of the latter equation into (9) we get the relaxed form of Eq. (3):

$$\begin{aligned} -\int _{\Omega }\rho \,c_V\,{\dot{T}}\,u\,\textrm{d}\Omega +\int _{\Omega }q\,u'\,\textrm{d}\Omega ={\tilde{q}}_{\ell }\,u(\ell )\,.\end{aligned}$$
(11)

Putting the weak form (8) together with the relaxed form (11), we arrived at a two-field mixed variational formulation which reads as follows. Find the duet \(T(t,\cdot )\in H^1(\Omega )\) and \(q(t,\cdot )\in L^2(\Omega )\) as trial functions satisfying a priori the BC (4) and the initial condition (6) such that

$$\begin{aligned} -\int _{\Omega }\rho \,c_V\,{\dot{T}}\,u\,\textrm{d}\Omega +\int _{\Omega }q\,u'\,\textrm{d}\Omega&={\tilde{q}}_{\ell }\,u(\ell )\quad \forall u\in H^1(\Omega )\,,\end{aligned}$$
(12)
$$\begin{aligned} \int _{\Omega }\frac{\tau }{\lambda }\,{\dot{q}}\,v\,\textrm{d}\Omega +\int _{\Omega }\frac{q\,v}{\lambda }\,\textrm{d}\Omega +\int _{\Omega }T'\,v\,\textrm{d}\Omega&=0\quad \forall v\in L^2(\Omega ) \end{aligned}$$
(13)

satisfying a priori the homogeneous form of the BC (4). Here, \(H^1(\Omega )\) stands for the Sobolev space of order 1 [39] and represents the space 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, as the first step of the derivation repeating the same process as presented in the latter subsection, i.e., relaxing Eq. (3) with the time-independent test function u(x), we obtain again the variational equation

$$\begin{aligned} -\int _{\Omega }\rho \,c_V\,{\dot{T}}\,u\,\textrm{d}\Omega +\int _{\Omega }q\,u'\,\textrm{d}\Omega ={\tilde{q}}_{\ell }\,u(\ell )\,.\end{aligned}$$
(14)

In the next step, after having divided the constitutive Eq. (2) by \(\lambda \), then multiplying it with the time-independent test function v(x) and integrating the product over the region \(\Omega \), again we have

$$\begin{aligned} \int _{\Omega }\left( \frac{\tau }{\lambda }\,{\dot{q}}+\frac{q}{\lambda }+\,T'-\frac{\kappa ^2}{\lambda }q''\,\right) v\,\textrm{d}\Omega =0\,, \end{aligned}$$
(15)

but now, integrating by parts the last term of this variational integral, we get

$$\begin{aligned}{} & {} -\int _{\Omega }\frac{\kappa ^2}{\lambda }q''\,v\,\textrm{d}\Omega =\\{} & {} \quad -\frac{\kappa ^2}{\lambda }\left[ q'\,v\right] _0^{\ell }+\int _{\Omega }\frac{\kappa ^2}{\lambda }q'\,v'\,\textrm{d}\Omega =\frac{\kappa ^2}{\lambda }\left[ q'(0)\,v(0)-q'(\ell )\,v(\ell )\right] +\int _{\Omega }\frac{\kappa ^2}{\lambda }q'\,v'\,\textrm{d}\Omega \,, \end{aligned}$$

the substitution of which into Eq. (15) results in the variational equation

$$\begin{aligned} \int _{\Omega }\frac{\tau }{\lambda }\,{\dot{q}}\,v\,\textrm{d}\Omega +\int _{\Omega }\frac{q\,v}{\lambda }\,\textrm{d}\Omega +\int _{\Omega }\,T'\,v\textrm{d}\Omega +\int _{\Omega }\frac{\kappa ^2}{\lambda }\,q'\,v'\,\textrm{d}\Omega =\frac{\kappa ^2}{\lambda }\left[ q'(\ell )\,v(\ell )-q'(0)\,v(0)\right] \,. \end{aligned}$$
(16)

As a result of this derivation, collecting Eqs. (14) and (16), again we arrive at a two-field mixed variational formulation but with a slightly different space regularity property from Eq. (13). Accordingly, we seek the variable pair \([T(t,\cdot ),\,q(t,\cdot )]\in H^1(\Omega )\) satisfying a priori the BC (4) and the initial condition (6) in such a way that

$$\begin{aligned}&-\int _{\Omega }\rho \,c_V\,{\dot{T}}\,u\,\textrm{d}\Omega +\int _{\Omega }q\,u'\,\textrm{d}\Omega ={\tilde{q}}_{\ell }\,u(\ell )\quad \forall u\in H^1(\Omega )\,,\end{aligned}$$
(17)
$$\begin{aligned}&\quad \int _{\Omega }\frac{\tau }{\lambda }\,{\dot{q}}\,v\,\textrm{d}\Omega +\int _{\Omega }\frac{q\,v}{\lambda }\,\textrm{d}\Omega +\int _{\Omega }\,T'\,v\,\textrm{d}\Omega +\int _{\Omega }\frac{\kappa ^2}{\lambda }\,q'\,v'\,\textrm{d}\Omega =\frac{\kappa ^2}{\lambda }\,q'(\ell )\,v(\ell )\nonumber \\&\quad -\frac{\kappa ^2}{\lambda }\,q'(0)\,v(0)\quad \forall v\in H^1(\Omega ) \end{aligned}$$
(18)

hold true, ensuring a priori the homogeneous form of the BC (4) once again.

2.2 hp-finite element discretizations

Let us now consider the hp-FE discretization of the domain \(\Omega \). Thenceforward, \(\Omega \) is divided into n physical sub-domain. Then, the master element \(\Omega _{\textrm{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 _{\textrm{h}}\) with the nodal points \(x_i\) and \(x_{i+1}\) of the FE mesh \(\Omega _{\textrm{h}}:0<x_1<x_2<\ldots<x_i<x_{i+1}<\ldots<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 piecewise 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\ge 2\), where \(L_k(\eta )\) are the orthogonal Legendre polynomials, \(k=3,4,\ldots ,p+1\) [35, 40, 41]. The independent trial-and-test functions, \(T,\;u\) and \(q,\;v\) are approximated on the i-th element of \(\Omega _{\textrm{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 space regularity assumptions on the trial functions T and q, as well as the test functions u and v in Eqs. (12, 13) for the MCV model and in Eqs. (17, 18) for the GK model. Here, p is the actual polynomial degree set to each element.

Table 1 Polynomial approximation spaces
Full size table

Let us represent now the variational Eqs (12, 13) and (17, 18) in matrix form. For the sake of simplicity, let \(^{\textrm{h}}{\varvec{t}}\) and \(^{\textrm{h}}{\varvec{q}}\) indicate the column matrices correspond to the time-dependent unknown temperature and heat flow coefficients, 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{aligned} \begin{array}{rl} {\mathcal {C}}\,^{\textrm{h}}\dot{{\varvec{t}}}+{\mathcal {Q}}^T\,^{\textrm{h}}{\varvec{q}}&{}={\varvec{d}}\,,\\ {\mathcal {T}}\,^{\textrm{h}}\dot{{\varvec{q}}}+{\mathcal {Q}}\,^{\textrm{h}}{\varvec{t}}+{\mathcal {K}}\,^{\textrm{h}}{\varvec{q}}&{}={\varvec{0}} \end{array} \end{aligned}$$
(19)

as semi-discrete equations which can be written in the simplified form

$$\begin{aligned} {\mathcal {A}}\,^{\textrm{h}}\dot{\varvec{\alpha }}+{\mathcal {B}}\,^{\textrm{h}}\varvec{\alpha }={\varvec{f}}\,, \end{aligned}$$
(20)

where

$$\begin{aligned} ^{\textrm{h}}\varvec{\alpha }=\left[ \begin{array}{c} ^{\textrm{h}}{\varvec{t}}\\ ^{\textrm{h}}{\varvec{q}} \end{array}\right] \,,\quad {\mathcal {A}}=\left[ \begin{array}{cc} {\mathcal {C}}&{}{\varvec{0}}\\ {\varvec{0}}&{}{\mathcal {T}} \end{array}\right] ,\quad {\mathcal {B}}=\left[ \begin{array}{cc} {\varvec{0}}&{}{\mathcal {Q}}^T\\ {\mathcal {Q}}&{}{\mathcal {K}} \end{array}\right] \quad \textrm{and}\quad {\varvec{f}}=\left[ \begin{array}{c} {\varvec{d}}\\ {\varvec{0}} \end{array}\right] , \end{aligned}$$
(21)

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.

It is important to note here that the time-dependence of the system is carried by the unknown coefficient values, i.e., the coefficients appeared in \(^{\textrm{h}}\varvec{\alpha }\) (\(\,^{\textrm{h}}{\varvec{t}}\) and \(\,^{\textrm{h}}{\varvec{q}}\)). In particular, the shape functions are not time-dependent at all, thereby seeking the numerical solution using the separation of variables technique. Even in the case of non-separable analytic solutions (if they exist), the solution can be approximated separately [15].

2.3 Time integration scheme for the semi-discrete system

The considered initial-boundary value problem as semi-discrete system contains seeking the function \(\varvec{\alpha }(t)\) (\(\,^{\textrm{h}}{\varvec{t}}(t)\) and \(\,^{\textrm{h}}{\varvec{q}}(t)\)) satisfying first-order time-dependent system (20) and the initial conditions (6) built-in \(\varvec{\alpha }_0\), for the numerical solution of which we use the generalized trapezoidal method. Accordingly, the time discretized version of Eq. (20) can be written as

$$\begin{aligned} ({\mathcal {A}}+\theta \,\Delta t\,{\mathcal {B}})\varvec{\alpha }^{t+\Delta t}=\left[ {\mathcal {A}}-(1-\theta )\Delta t\,{\mathcal {B}}\right] \varvec{\alpha }^t+\Delta t\left[ \theta {\varvec{f}}^{t+\Delta t}+(1-\theta ){\varvec{f}}^t\right] \,,\end{aligned}$$
(22)

which can be simplified into the form

$$\begin{aligned} {\mathcal {E}}\varvec{\alpha }^{t+\Delta t}={\mathcal {D}}\varvec{\alpha }^t+\Delta t\left[ \theta {\varvec{f}}^{t+\Delta t}+(1-\theta ){\varvec{f}}^t\right] \,,\end{aligned}$$
(23)

where \(\Delta t\) is the time step size assumed to be constant and the parameter \(\theta \) is taken to be in the interval [0, 1], as well as \({\mathcal {E}}={\mathcal {A}}+\theta \,\Delta t\,{\mathcal {B}}\) and \({\mathcal {D}}={\mathcal {A}}-(1-\theta )\Delta t\,{\mathcal {B}}\). Thus, the amplifier matrix of the resulting system

$$\begin{aligned} \varvec{\alpha }^{t+\Delta t}=\mathcal {{\tilde{A}}}(\theta )\,\varvec{\alpha }^t+\varvec{{\tilde{f}}}(\theta ) \end{aligned}$$
(24)

is computed as \(\varvec{{\tilde{A}}}={\mathcal {E}}^{-1}{\mathcal {D}}\) (\(\varvec{{\tilde{f}}}=\Delta t[\theta \,{\mathcal {E}}^{-1}{\varvec{f}}^{t+\Delta t}+(1-\theta ){\mathcal {E}}^{-1}{\varvec{f}}^t]\)). During the numerical investigations, we choose 1 for the parameter \(\theta \) in order to make the algorithm unconditionally stable. Besides, on each convergence test and simulation, the spectral radius and the eigenvalues of \(\mathcal {{\tilde{A}}}\), i.e., the stability requirements on \(\mathcal {{\tilde{A}}}\) [13, 15] is computationally checked.

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 analysis. 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 is measured, respectively, in the maximum norm of the temperature and the heat flow

$$\begin{aligned} e^{\textrm{hp}}=\frac{\underset{t\in [t_0,t_1]}{\max }\left| \,^{\textrm{h}}{\varvec{t}}-{\varvec{t}}_{\textrm{ref}}\right| }{\underset{t\in [t_0,t_1]}{\max }\left| {\varvec{t}}_{\textrm{ref}}\right| }\quad \text {and}\quad e^{\textrm{hp}}=\frac{\underset{t\in [t_0,t_1]}{\max }\left| \,^{\textrm{h}}{\varvec{q}}-{\varvec{q}}_{\textrm{ref}}\right| }{\underset{t\in [t_0,t_1]}{\max }\left| {\varvec{q}}_{\textrm{ref}}\right| }. \end{aligned}$$
(25)

The related thermodynamic system (20) is solved numerically for the relatively small time interval \(t\in [0,10]\) s, using the implicit time integration scheme, i.e., the setting \(\theta =1\) [13, 15]. The number of the constant time steps 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 [9] to be realistic; therefore, we choose \(\tau =\{0.05, 0.15, 0.3\}\) s, with \(\kappa ^2=8\cdot 10^{-6}\) m\(^2\).

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

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

$$\begin{aligned} {\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 \textrm{and}\quad {\tilde{q}}_{\ell }(t)=0\end{aligned}$$
(26)

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 an equilibrium state. The \(c_1\) and \(c_2\) coefficients are chosen based on [4] to keep the excitation as realistic as possible.

During the p-convergence tests, 8-, 52- and 20-element equidistant meshes are fixed on the domain \(\Omega \) for the GK model along with the parameter settings \(\kappa ^2=0.8\) m\(^2\), \(\kappa ^2=0.000008\) m\(^2\) and the MCV model, respectively, while the polynomial degree p is ranging from 2 to 8 for all of FEs. 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\) m\(^2\), \(\kappa ^2=0.000008\) m\(^2\) and the MCV model. In each mesh refinement step, the polynomial degree p is equal to 2 and unchanged for all elements. Furthermore, the number of degrees of freedom (DOF) is equal to the number of the unknown coefficients occurring in \(\varvec{\alpha }\), i.e., the size of \(\varvec{\alpha }\).

Fig. 2
figure 2

Convergence histories of the relative errors measured in the maximum norm for the front- and rear-side temperature, as well as the mid-side heat flow—MCV model

Full size image
Fig. 3
figure 3

Convergence histories of the relative errors measured in the maximum norm for the front- and rear-side temperature, as well as the mid-side heat flow—GK model with \(\kappa ^2=0.8\;\textrm{m}^2\)

Full size image
Fig. 4
figure 4

Convergence histories of the relative errors measured in the maximum norm for the front- and rear-side temperature, as well as the mid-side heat flow—GK model with \(\kappa ^2=0.000008\;\textrm{m}^2\)

Full size image

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 m\(^2\) and 0.000008 m\(^2\)) are plotted against the number of DOF on log–log scales in Figs. 2 and 3, 4 for the relaxation time \(\tau =0.3\) s, \(\tau =0.15\) s, and \(\tau =0.05\) s, respectively.

The relative error behavior exhibits fast convergences not only for p- but also for h-approximation. 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 a 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 an algebraic type of convergence instead.

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 \) decreases, achieving a lower relative error value. Besides, it can be seen that the error level is around a very small value (approximately \(10^{-8}\)), i.e., the computational accuracy is very high in Figure 6.

Besides, it is worth mentioning that during all of the convergence tests, the amplifier matrix has been spectrally stable as expected in Figure 10.

Fig. 5
figure 5

Time series of the dimensionless rear-side temperature for the MCV and GK model with \(\kappa ^2=0.000008\;\textrm{m}^2\) – \(\tau =0.3\) s, \(p=10\) and \(n=100\)

Full size image
Fig. 6
figure 6

Time series of the dimensionless rear-side temperature for the MCV and GK model with \(\kappa ^2=0.000008\;\textrm{m}^2\) – \(\tau =0.15\) s, \(p=10\) and \(n=100\)

Full size image
Fig. 7
figure 7

Time series of the dimensionless rear-side temperature for the MCV and GK model with \(\kappa ^2=0.000008\;\textrm{m}^2\) – \(\tau =0.05\) s, \(p=10\) and \(n=100\)

Full size image
Fig. 8
figure 8

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\;\textrm{m}^2\)—\(p=10\) and \(n=100\)

Full size image
Fig. 9
figure 9

Time series of the dimensionless front-side temperature for the MCV and GK model with \(\kappa ^2=0.000008\;\textrm{m}^2\) – \(\tau =0.3\) s, \(p=10\) and \(n=100\)

Full size image
Fig. 10
figure 10

Time series of the dimensionless front-side temperature for the MCV and GK model with \(\kappa ^2=0.000008\;\textrm{m}^2\) – \(\tau =0.15\) s, \(p=10\) and \(n=100\)

Full size image
Fig. 11
figure 11

Time series of the dimensionless front-side temperature for the MCV and GK model with \(\kappa ^2=0.000008\;\textrm{m}^2\) – \(\tau =0.05\) s, \(p=10\) and \(n=100\)

Full size image
Fig. 12
figure 12

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\;\textrm{m}^2\) – \(p=10\) and \(n=100\)

Full size image

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 m\(^2\)), the time histories of the dimensionless rear- and front-side temperature are depicted, separately for the three different relaxation times \(\tau =0.3\) s, \(\tau =0.15\) s, and \(\tau =0.05\) s, in Figs. 5, 7 and, 9, 11, 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. Figures 8 and 12 represent the time series of the dimensionless rear- and front-side temperature for “over-diffuse” thermodynamic system, i.e., for \(\kappa ^2=0.8\) m\(^2\), zooming the hp-FE solutions in the rapidly change and very small time region [0, 0.02] s. These figures show that the \(\tau \)-value has no significant effect on the solution for “over-diffuse” system.

It is important to mention here that oscillation-free simulation results of the MCV model are regained from the GK model as well if we precisely substitute zero into the GK model for the diffusivity parameter \(\kappa ^2\). From both the GK and the MCV model, the simulation results of the Fourier model can be obtained without oscillations, substituting precisely zero values for both the relaxation time \(\tau \) and the diffusivity parameter \(\kappa ^2\) in the GK and/or the MCV model. Namely, the hp-FEs exhibit oscillation-free behavior even if \(\tau \) and \(\kappa ^2\) are exactly equal to zero.

4 Discussion

Here, we have presented novel numerical approaches which are based on (1) newly derived two-field variational formulations that consider the temperature and heat flux field as independent variables and (2) new, advanced hp-version mixed FEMs. Within this new concept, we have focused specifically on the FE solution of the Guyer–Krumhansl equation, as that model could have a 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, contrary to COMSOL [11]. That is essential since the usual approaches do not work, and it is not possible to use the temperature as a single field variable for heat flux BCs [2]. The advantage of the presented hp-FEMs is as follows:

  1. (i)

    Since the derived variational formulations and the related FEMs are two-field ones, we obtain numerical results directly for both the temperature field and the heat flux, thereby decreasing the computational cost and time of the postprocessing subroutine of the code.

  2. (ii)

    Fast p-convergence with high precision has experienced in both the temperature and the heat flow computations. It means that highly accurate FE solutions can be achieved with relatively low-density mesh by means of very high degree polynomial approximation, i.e., the desired (high) accuracy is mainly reached much earlier by the p-type approximation because the p-convergence is much faster than the h-convergence, thereby verifying the effectiveness of the hp-FE code.

  3. (iii)

    The developed hp-FEMs are locking-free, producing robust (uniformly stable) FE solutions at each time instant, i.e., providing reliable numerical results for both the temperature and the heat flux.

The only disadvantage of the constructed mixed hp-FEMs is that these can lead to large sparse matrices, thereby increasing (but not significantly) the memory usage and the runtime during the solution. However, this can be easily remedied by the generation and use of preconditioners for the main matrix blocks. However, the preconditioning has not been necessary up to now, thanks to the excellent convergence test results described above.

Moreover, that technique is incredibly more efficient than the one offered by COMSOL and runs about four 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 massive increase in speed would enable the real-time monitoring of heterogeneous materials with complex inner structures without massive computational capacity.