New perspectives for modelling ballistic-diffusive heat conduction

Abstract

Among the three heat conduction modes (diffusive, second sound, and ballistic propagation), the ballistic one is the most difficult to model. In the present paper, we discuss its physical interpretations and show different alternatives to its modelling. We highlight two of them: a thermo-mechanical model in which the generalized heat equation—the so-called Maxwell–Cattaneo–Vernotte equation—is coupled to thermal expansion. The other one uses internal variables in the framework of non-equilibrium thermodynamics. For the first one, we developed a numerical solution and tested it on the heat pulse experiment performed by McNelly et al. on NaF samples. For the second one, we found an analytical solution that emphasizes the role of boundary conditions. This analytical method is used to validate an earlier developed numerical code.

Appendices

Appendix A. Stability of the numerical scheme

Since the scheme that we utilize here for demonstration is at least partially explicit, we must investigate the corresponding stability properties. Following [12, 56,57,58], we use the von Neumann method and the Jury conditions, which is suitable for linear equations. However, due to the appearance of the nonlinear term $(T-T_0)\partial _xv$, we have to do first the following. Estimating the maximum of the temperature field T is making the nonlinear term to be linear with substituting a constant for T [38]. Unfortunately, it is an open mathematical question of how to make the proper estimation for T. Here, we assume that $C=\max _{(j,n)} (T^n_j-1)$, for all x and t is 1. Next, what serves for the starting point is the discrete version of the dispersion relation, i.e. assuming the solution for the difference equations in the form of a plane wave:

$$\begin{aligned} u_j^n=\xi ^n e^{ikj\Delta x}, \end{aligned}$$

(25)

in which $\xi$ is called growth factor; this is the amplitude of the wave and must be bounded from above for stability: $|\xi |\le 1$. Furthermore, i is the imaginary unit, and k is the wavenumber. After substituting Eq. (25) into (15), we obtain the system of linear algebraic equations:

$$\begin{aligned} {\mathbf {M}}\cdot \left( \begin{array}{llll} T_0&\varepsilon _0&q_0&v_0 \end{array}\right) ^T =0, \end{aligned}$$

(26)

where the coefficient matrix is

$$\begin{aligned} {\mathbf {M}}= \left( \begin{array}{llll} \xi -1&{} 0 &{} \frac{\Delta t}{\Delta x} P_1\left( e^{ik\Delta x}-1\right) &{} \frac{\Delta t}{\Delta x} P_2 C \left( e^{ik\Delta x}-1\right) \\ 0 &{} \xi -1 &{}0 &{} - \frac{\Delta t}{\Delta x} \left( e^{ik\Delta x}-1\right) \\ \frac{\Delta t}{\Delta x} \frac{\lambda _d}{\tau _d} \xi \left( 1-e^{-ik\Delta x}\right) &{}0 &{}\xi -1+\frac{\Delta t}{\Delta x} &{}0\\ \frac{\Delta t}{\Delta x} \alpha \xi \left( 1-e^{-i k\Delta x}\right) &{}- \frac{\Delta t}{\Delta x} \xi \left( 1-e^{-ik\Delta x}\right) &{}0 &{}\xi -1 \end{array} \right) . \end{aligned}$$

Then the characteristic polynomial for $\xi$ can be formally expressed as

$$\begin{aligned} p(\xi )=a_4\xi ^4+a_3\xi ^3+a_2\xi ^2+a_1\xi +a_0, \end{aligned}$$

(27)

where the coefficients are found to be

$$\begin{aligned} a_4&=1; \nonumber \\ a_3&=\frac{\Delta t \Delta x^2 - 2(\cos (k\Delta x)-1)\Delta t (P_1 \lambda _d + \tau _d + P_2 C \alpha \tau _d)\Big )}{\Delta x^2 \tau _d}-4; \nonumber \\ a_2&=\frac{4(\cos (k\Delta x)-1)^2 \Delta t^4 P_1 \lambda _d -3 \Delta x^4 (\Delta t - 2 \tau _d)}{\Delta x^4\tau _d}-\nonumber \\&-\frac{2(\cos (k\Delta x)-1)\Delta t^2 (\Delta t + \Delta t P_2 C \alpha -2 P_1 \lambda _d -2 \tau _d - 2 P_2 C \alpha \tau _d)}{\Delta x^2\tau _d};\nonumber \\ a_1&=\frac{(3\Delta t - 4 \tau _d)\Delta x^2 +2(\cos (k\Delta x)-1)\Delta t^2 \Big (\Delta t + P_2 C \alpha - P_1 \lambda _d -(1+ P_2 C \alpha \tau _d )\Big )}{\Delta x^2 \tau _d} \nonumber \\ a_0&=1-\Delta t/\tau _d. \end{aligned}$$

(28)

The roots of that characteristic polynomial can be restricted into the unit circle on a complex plane, using Jury conditions [58], i.e.

• $p(\xi =1) \ge 0$;

• $p(\xi =-1) \ge 0$;

• $|a_0|\le |a_4|$;

• $|b_0 |>|b_3|$, where $b_k= \left| \begin{array}{cc} a_0 &{} a_{n-k} \\ a_n &{} a_k \end{array} \right|$, with $n =4$ and $k=3$;

• $|c_0 |>|c_2|$, where $c_k= \left| \begin{array}{cc} b_0 &{} b_{n-1-k} \\ b_{n-1} &{} b_k \end{array} \right|$, with $n =4$ and $k= 2$.

Unfortunately, the analytical form of these inequalities does not allow to simplify the expressions reasonably; therefore, we decided to check them numerically. However, let us remark that the stability conditions are in agreement with the thermodynamic restrictions, including the relations between the material parameters.

Appendix B. Analytical solutions

The complete solution of Eq. (23) is split into two parts. The first one lasts until the end of the heat pulse $\tau _\Delta$, while the next one is having a constant boundary, independently of the time.

1.1 Section I ($0<t\le t_\Delta$)

Starting with the first section, let us assume that we can write the heat flux in the form $q=v+w$, in which w serves to separate the time-dependent inhomogeneous boundary from v which becomes the homogeneous solution of Eq. (23), and

$$\begin{aligned} q=v+w; \quad w = \Big ( 1-\frac{x}{L} \Big ) q_0(t). \end{aligned}$$

(29)

Substituting Eq. (29) into (23), and executing the derivatives, we obtain the following inhomogeneous partial differential equation:

$$\begin{aligned} a \dddot{v}+b\ddot{v}+{\dot{v}}=v''+c{\dot{v}}''-f(x,t), \end{aligned}$$

(30)

where $f(x,t)=a\dddot{w}+b\ddot{w}+{\dot{w}}$. Using the splitting $v=q-w$, the initial conditions for v are not completely equal to zero, that is,

$$\begin{aligned} v(x,t=0)={\dot{v}}(x,t=0)=0, \quad \ddot{v}(x,t=0) = \frac{4\pi ^2 x}{\tau _{\Delta }^2L}. \end{aligned}$$

(31)

Applying the classical procedure of separation of variables, we seek the solution in the form $v(x,t)=\phi (t)X(x)$, in which $\phi (t)$ describes the time evolution, and X(x) stands for the spatial behaviour. Furthermore, following [64], we assume that the inhomogeneity f(x, t) can be represented in the space spanned by the eigenfunctions of Laplacian operator. Thus, we determine the eigenspace for the homogeneous part ($f(x,t)=0$),

$$\begin{aligned} \frac{a\dddot{\phi }+b\ddot{\phi }+{\dot{\phi }}}{\phi +c{\dot{\phi }}}=\frac{X''}{X}, \end{aligned}$$

(32)

where $X''/X = \beta$ yields the eigenvalues and the eigenfunctions, with the homogeneous boundaries $X(0)=X(L)=0$:

$$\begin{aligned} X_n(x)=\sin (\sqrt{\beta _n}x), \end{aligned}$$

(33)

with $\beta _n = \Big ( \frac{n\pi }{L}\Big )^2$ being the eigenvalues. Then, exploiting the linearity of Eq. (23), v(x, t) reads

$$\begin{aligned} v(x,t)=\sum _{n=1}^{\infty } \phi _n(t)\sin \Big (\frac{n\pi }{L}x\Big ). \end{aligned}$$

(34)

Now, we can substitute Eq. (34) into (30):

$$\begin{aligned} \sum _{n=1}^{\infty }\Big [ a\dddot{\phi }_n + b \ddot{\phi }_n + (1+c\beta _n){\dot{\phi }}_n + \beta _n \phi _n \Big ]\sin \Big ( n\pi \frac{x}{L}\Big ) = -f(x,t), \end{aligned}$$

(35)

which equation requires the Fourier sin series expansion of f(x, t), i.e. $f(x,t)=\sum _{n=1}^{\infty }f_n(t)\sin \Big ( n\pi \frac{x}{L}\Big )$, and $f_n(t)$ is expressed using Eqs. (24) and (29):

$$\begin{aligned}&{\displaystyle f_n(t)=\frac{2}{L}\Bigg [\Big (\frac{2\pi }{\tau _{\Delta }} -a \frac{8\pi ^3}{\tau _{\Delta }^3}\Big ) \sin \Big ( 2\pi \frac{t}{\tau _{\Delta }}\Big ) + b\frac{4\pi ^2}{\tau _{\Delta }^2}\cos \Big ( 2\pi \frac{t}{\tau _{\Delta }}\Big ) \Bigg ] }\\ \nonumber&\quad {\displaystyle \times \int _{0}^{L} \Big ( 1-\frac{x}{L}\Big ) \sin \Big ( \frac{n\pi }{L}x\Big ) \text {d}x. } \end{aligned}$$

(36)

After integration and simplification, we obtain the nth term as follows:

$$\begin{aligned} f_n(t) = \frac{2}{n\pi }\Bigg [ \Big ( \frac{2\pi }{\tau _{\Delta }}-\frac{8a\pi ^3}{\tau _{\Delta }^3}\Big ) \sin \Big ( 2\pi \frac{t}{\tau _{\Delta }} \Big ) + \frac{4b\pi ^2}{\tau _{\Delta }^2}\cos \Big ( 2\pi \frac{t}{\tau _{\Delta }}\Big ) \Bigg ]. \end{aligned}$$

(37)

Now, using Eq. (37), Eq. (35) shall lead to the solution of the time evolution for v, i.e. $\phi _n(t)$ is prescribed as

$$\begin{aligned} a \dddot{\phi }_n +b\ddot{\phi }_n+(1+c\beta _n){\dot{\phi }}_n+\beta _n\phi _n=-f_n(t), \end{aligned}$$

(38)

and its general solution is

$$\begin{aligned} \phi _n = Y_{1_n} e^{x_1 t} + Y_{2_n} e^{x_2 t} + Y_{3_n} e^{x_3 t} + A \sin (g t) + B \cos (g t), \end{aligned}$$

(39)

where $x_{1,2,3}$ are the roots of the characteristic polynomial $ax^3+bx^2+(1+c\beta )x+\beta =0$, $g=2\pi / \tau _{\Delta }$, and $A, \ B$ are the parameters of the particular solution $\phi _p$. Now substituting $\phi _p=A\sin (gt)+B\cos (gt)$ into Eq. (38),

$$\begin{aligned} a \dddot{\phi }_p + b \ddot{\phi }_p + (1+c\beta _n){\dot{\phi }}_p+\beta _n\phi _p = K_1 \sin (gt) + K_2\cos (gt), \end{aligned}$$

(40)

with

$$\begin{aligned} K_1 = -\frac{2}{n \pi }\Bigg ( \frac{2\pi }{\tau _{\Delta }} - \frac{8a\pi ^3}{\tau _{\Delta }^3} \Bigg ); \quad K_2 = - \frac{2}{n\pi }\frac{4b\pi ^2}{\tau _{\Delta }^2}, \end{aligned}$$

(41)

as well as

$$\begin{aligned} -aAg^3-bBg^2+A(1+\beta _nc)g+\beta _n B&=K_2, \nonumber \\ aBg^3-bAg^2-Bg(1+\beta _nc)+\beta _n A&=K_1, \end{aligned}$$

(42)

which system gives the expressions for A and B. In order to determine the coefficients $Y_{1_n,2_n,3_n}$, we must take the initial conditions into account:

$$\begin{aligned} \phi (0)&= Y_{1_n}+Y_{2_n}+Y_{3_n}+B=0, \nonumber \\ {\dot{\phi }}(0)&= x_1Y_{1_n}+x_2Y_{2_n}+x_3Y_{3_n}+Ag=0, \nonumber \\ \ddot{\phi }(0)&= x_1^2Y_{1_n}+x_2^2Y_{2_n}+x_3^2Y_{3_n}-Bg^2=\frac{4\pi ^2 x}{\tau _{\Delta }^2L}. \end{aligned}$$

(43)

Formally, we obtained the solution of the heat flux q for the first session that lasts until $\tau _\Delta$,

$$\begin{aligned} q_I(x,t)= & {} \sum _{n=1}^{\infty }\Big ( Y_{1_n} e^{x_1 t} + Y_{2_n} e^{x_2 t} + Y_{3_n} e^{x_3 t} + A \sin (g t) + B \cos (g t) \Big ) \sin (\frac{n\pi }{L}x) \nonumber \\&+ \Big (1-\frac{x}{L} \Big ) q_0(t) . \end{aligned}$$

(44)

We can use this expression to reconstruct the temperature field, utilizing the balance of internal energy:

$$\begin{aligned}&{\displaystyle T_I(x,t) =\frac{t}{L\tau _{\Delta }} - \frac{1}{2\pi L}\sin (gt) } \nonumber \\&\quad {\displaystyle - \sum _{n=1}^{\infty }\frac{n\pi \cos (\frac{n\pi }{L}x)}{L \tau _{\Delta }} \Bigg ( \Bigg [ \sum _{j=1}^3 \frac{e^{x_jt}-1}{x_j}Y_{jn} \Bigg ] + \frac{B_n}{g}\sin (gt)+\frac{A_n-A_n\cos (gt)}{g}\Bigg )}. \end{aligned}$$

(45)

1.2 Section II ($t \ge \tau _\Delta$)

Now turning our attention to the second session of the solution (in which $t\ge \tau _\Delta$), we introduce a new variable for time such as ${\tilde{t}}=t-\tau _{\Delta }$ for convenience. Thus, the initial conditions become

$$\begin{aligned} q_{II}(x,{\tilde{t}}=0)&= q_I(x,t=\tau _{\Delta }), \nonumber \\ {\dot{q}}_{II}(x,{\tilde{t}}=0)&= {\dot{q}}_I(x,t=\tau _{\Delta }), \nonumber \\ \ddot{q}_{II}(x,{\tilde{t}}=0)&= \ddot{q}_I(x,t=\tau _{\Delta }). \end{aligned}$$

(46)

Since the inhomogeneity of f(x, t) disappears as $w\equiv 0$, we can immediately begin with

$$\begin{aligned} a \dddot{\gamma }_n +b\ddot{\gamma }_n+(1+c\beta _n){\dot{\gamma }}_n+\beta _n\gamma _n=0. \end{aligned}$$

(47)

The roots of the characteristic polynomial $x_{1,2,3}$ remain unchanged, and the general solution reads

$$\begin{aligned} \gamma _n(t)=C_{1_n}e^{x_1t}+C_{2_n}e^{x_2t}+C_{3_n}e^{x_3t}, \end{aligned}$$

(48)

in which the unknown $C_{1,2,3}$ coefficients are determined using the initial conditions (46):

$$\begin{aligned} C_{1_n}+C_{2_n}+C_{3_n}=\phi _n(t=\tau _{\Delta }), \nonumber \\ x_1C_{1_n}+x_2C_{2_n}+x_3C_{3_n}={\dot{\phi }}_n(t=\tau _{\Delta }), \nonumber \\ x_1^2C_{1_n}+x_2^2C_{2_n}+x_3^2C_{3_n}=\ddot{\phi }_n(t=\tau _{\Delta }). \end{aligned}$$

(49)

Again, the temperature distribution is recovered using the balance equation.

$$\begin{aligned} T_{II}(x,{\tilde{t}}) = - \sum _{n=1}^{\infty }\frac{n\pi \cos (\frac{n\pi }{L}x)}{L \tau _{\Delta }} \Bigg [ \sum _{j=1}^3 \frac{e^{x_j{\tilde{t}}}-1}{x_j}C_{j_n} \Bigg ] + T_{I}(x,t=\tau _{\Delta } ) , \end{aligned}$$

(50)

where ${\tilde{t}}=t-\tau _{\Delta }$, and $T_{I}(x,t=\tau _{\Delta })$ is defined by Eq. (45).

1.3 Calculating Q(x, t)

The heat current density Q(x, t) can be calculated from Eq. (18) using the previously determined T(x, t) and q(x, t) functions. The solution again is divided into two parts, namely $Q_{I}(x,t)$ and $Q_{II}(x,t)$, where the $T_{I,II}(x,t)$ and $q_{I,II}(x,t)$ functions are used accordingly. The necessary derivatives needed for Q(x, t) are the following:

$$\begin{aligned} \partial _t q_{I}(x,{\hat{t}})= & {} \sum _{n=1}^{\infty } \Bigg [ \sum _{k=1}^{3}Y_{kn}x_{k}e^{x_k t} + A_ng\cos (gt) - B_ng\sin (gt)\Bigg ]\sin \Big ( \frac{n\pi }{L}x\Big ) \nonumber \\&{\displaystyle + \Big ( 1-\frac{x}{L} \Big ) g \sin (gt)}, \end{aligned}$$

(51)

$$\begin{aligned} \displaystyle \partial _x T_{I}(x,{\hat{t}})= & {} \sum _{n=1}^{\infty } \frac{(n\pi )^2\sin \Big ( \frac{n\pi }{L}x \Big )}{L^2\tau _{\Delta }}\Bigg [ \sum _{k=1}^3 \frac{e^{x_k t}-1}{x_k}Y_{kn} + \frac{B_n}{g}\sin (gt) \nonumber \\&{\displaystyle + \frac{A_n-A_n\cos (gt)}{g} \Bigg ]}, \end{aligned}$$

(52)

$$\begin{aligned} \partial _t q_{II}(x,{\hat{t}})= & {} \sum _{n=1}^{\infty } \sum _{k=1}^{3}C_{kn}x_{k}e^{x_k{\hat{t}}}\sin \Big ( \frac{n\pi }{L}x\Big ), \end{aligned}$$

(53)

$$\begin{aligned} \partial _x T_{II}(x,{\hat{t}})= & {} \sum _{n=1}^{\infty } \frac{(n\pi )^2\sin \Big ( \frac{n\pi }{L}x \Big )}{L^2\tau _{\Delta }}\Bigg [ \sum _{k=1}^3 \frac{e^{x_k {\hat{t}}}-1}{x_k}C_{kn} + \sum _{k=1}^3 \frac{e^{x_k \tau _{\Delta }}-1}{x_k}Y_{kn} \nonumber \\&+ \frac{B_n}{g}\sin (g\tau _{\Delta }) + \frac{A_n-A_n\cos (g\tau _{\Delta })}{g} \Bigg ]. \end{aligned}$$

(54)

Putting everything together the heat flux density for the first time interval is:

$$\begin{aligned} Q_{I}(x,t)=\frac{1}{\kappa } \sum _{n=1}^{\infty } \frac{L\cos \Big ( \frac{n\pi }{L}x\Big )}{n\pi }{\mathcal {F}}(A_n,B_n,t,x_k) -\Big ( x-\frac{x^2}{2L}\Big )\Big ( 1-\cos (gt)+\tau _qg\sin (gt)\Big ), \end{aligned}$$

(55)

where

$$\begin{aligned} {\mathcal {F}}(A_n,B_n,t,x_k)= & {} A_n\sin (gt)+B_n\cos (gt)+\sum _{k=1}^3 Y_{kn}e^{x_kt} + \tau _q A_ng\cos (gt)-\tau _qB_ng\sin (gt) \nonumber \\&+ \tau _q\sum _{k=1}^3Y_{kn}x_k e^{x_kt} + \frac{n^2\pi ^2 B_n}{L^2 g}\sin (gt) + \frac{n^2\pi ^2}{L^2}\frac{A_n-A_n \cos (gt)}{g} +\frac{n^2\pi ^2}{L^2}\sum _{k=1}^3 \frac{e^{x_kt}-1}{x_k}Y_{kn} . \end{aligned}$$

(56)

Similarly, the heat flux density for the second interval will be:

$$\begin{aligned} Q_{II}(x,t) = \frac{1}{\kappa }\sum _{n=1}^{\infty }\frac{L\cos \Big ( \frac{n\pi }{L}x\Big )}{n\pi }{\mathcal {G}}(A_n,B_n,x_k,{\hat{t}},\tau _{\Delta }), \end{aligned}$$

(57)

in which

$$\begin{aligned} {\mathcal {G}}(A_n,B_n,x_k,{\hat{t}},\tau _{\Delta })= & {} \sum _{k=1}^{3}C_{kn}e^{x_k{\hat{t}}} + \tau _q\sum _{k=1}^3 C_{kn}x_ke^{x_k{\hat{t}}}+\frac{n^2\pi ^2}{L^2}\Bigg [\sum _{k=1}^3 \frac{e^{x_k{\hat{t}}}-1}{x_k}C_{kn} \nonumber \\&+ \sum _{k=1}^3 \frac{e^{x_k\tau _{\Delta }}-1}{x_k}Y_{kn} + \frac{B_n}{g}\sin (g\tau _{\Delta })+ \frac{A_n-A_n\cos (g\tau _{\Delta })}{g} \Bigg ] . \end{aligned}$$

(58)