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

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.


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: in which q, λ and T are the heat flux, thermal conductivity and temperature field; Π describes the dynamic pressure as a consequences of the presence of velocity gradient ∂ x v; and the stress tensor σ is proportional with the deformation ε.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, in which the partial time derivative of the heat flux ( q) appears with a coefficient τ called relaxation time.Moreover, Eq. ( 2) consists of the second-order spatial derivative of q as well, denoted by q ′′ .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 κ 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 κ 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, where ρ 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 as Dirichlet-and Neumann-type BCs at x = 0 and x = ℓ, respectively, as well as the temporal prescriptions T (0, x) = T 0 (x) and q(0, x) = q 0 (x) (6) as initial conditions at t 0 = 0 s.

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].

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 Ω with setting κ 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 ∈ H 1 (Ω) and q ∈ L 2 (Ω) as trial functions satisfying a priori the BC (4) and the initial condition ( 6) such that satisfying a priori the homogeneous form of the BC (4).Here, H 1 (Ω) stands for the Sobolev space of order 1 [23] and represents the regularity property for T and u while L 2 (Ω) defines square integrable function space for q and v.

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) ∈ H 1 (Ω) satisfying a priori the BC (4) and the initial condition (6) in such a way that hold true, ensuring a priori the homogeneous form of the BC (4) once again.Notably, again, for both variables (T, q) ∈ H 1 (Ω), which stands as a crucial difference for other FE approaches for heat equations.

hp-finite element discretizations
Let us consider now the hp-FE discretization of the domain Ω.Thenceforward Ω is divided into n physical subdomain.Then the master element 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 where L k (η) are the orthogonal Legendre polynomials, k = 3, 4, . . ., p + 1 [18,24,25].The independent trial-and test functions, T, u and q, v are approximated on the i-th element of Ω 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.
Let us represent now the variational equations ( 7)-( 8) and ( 9)-( 10) in matrix form.For the sake of simplicity let h t and h q indicate the column matrices corresponding to the time-dependent unknown coefficients of the temperature and

Model
Function T, u q, v 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 which can be written in the simplified form where in which C, T and K denote the consistent matrices of the specific heat, the relaxation term and the heat conductivity, respectively, while Q is the consistent coupling matrix of the system.

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 The related thermodynamic system ( 12) is solved numerically for the relatively small time interval t ∈ [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 ∆t = 0.001 s.The considered bodies of length ℓ = 0.005 m are made of rock-like materials which have c V = 800 (J/kgK), ρ = 2600 kg/m 3 .Additionally, the new parameters τ and κ 2 are determined based on earlier experimental data [5] to be realistic, therefore we choose τ = {0.05,0.15, 0.3} s, with κ 2 = 8 • 10 −6 m 2 .However, in order to challenge the numerical procedure, we also test for unrealistically high κ 2 = 0.8 m 2 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 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 Ω for the GK model with the parameter settings κ 2 = 0.8 m 2 , κ 2 = 0.000008 m 2 and the MCV model, respectively, while the polynomial degree p is ranging from 2 to 8 for each element.During the hconvergence studies, the domain Ω 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 κ 2 = 0.8 m 2 , κ 2 = 0.000008 m 2 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 α.
The relative error-convergence curves of the front-and rear side temperature, T (t, 0) and T (t, ℓ), as well as the mid-side heat flow q(t, ℓ/2) obtained for the MCV model and the GK model with relatively large and small value of κ 2 (0.8 m 2 and 0.000008 m 2 ) are plotted against the number of DOF on log-log scales in Figs. 1 and 2-3 for the relaxation time τ = 0.3 s, τ = 0.15 s and τ = 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 higherorder 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 τ on the convergence rates.However, the convergence curves are shifted down (without changing their slopes) as τ 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 κ 2 (0.000008 m 2 ), the histories of the dimensionless rear-and front side temperatures are depicted, separately for three different relaxation times τ = 0.3 s, τ = 0.15 s and τ = 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 κ 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 exhibit very well that the τ-value has no significant effect on the solution for "over-diffuse" system.

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.

Figure 1 :Figure 4 :Figure 5 :Figure 6 :
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

Figure 7 :Figure 8 :Figure 9 :Figure 10 :
Figure 7: Influence of τ on the time series of the dimensionless rear side temperature for the GK model with the over-diffuse setting κ 2 = 0.8 m 2p = 10 and n = 100

Figure 11 :
Figure 11: Influence of τ on the time series of the dimensionless front side temperature for the GK model with the over-diffuse setting κ 2 = 0.8 m 2 p = 10 and n = 100