On the propagation of temperature-rate waves and traveling waves in rigid conductors of the Graffi--Franchi--Straughan type

TL;DR

This paper investigates second-sound phenomena in Graffi--Franchi--Straughan type conductors, revealing nonlinear effects like wave blow-up and unique waveforms due to temperature-dependent relaxation times, using analytical and numerical methods.

Contribution

It introduces a detailed analysis of temperature-rate and traveling waves in a new class of conductors with nonlinear effects caused by temperature-dependent relaxation times.

Findings

Finite-time temperature-rate wave blow-up observed.

Thermal traveling waves exhibit a 'tongue' waveform.

Nonlinear effects are significant in these conductors.

Abstract

We examine second-sound phenomena in a class of rigid, thermally conducting, solids that are described by a special case of the Maxwell--Catteneo flux law. Employing both analytical and numerical methods, we examine both temperature-rate waves and thermal traveling waves in this class of thermal conductor, which have recently been termed Graffi--Franchi--Straughan type conductors. In the present study, the temperature-dependent nature of the thermal relaxation time, which is the distinguishing feature of this class of conductors, gives rise to a variety of nonlinear effects; in particular, finite-time temperature-rate wave blow-up and temperature traveling waveforms which exhibit a "tongue". The presentation concludes with a discussion of possible follow-on studies.

Figures

Equations

$$\mathbf{q}+\alpha \vartheta \mathcal{K}(\vartheta ){\mathbf{q}}_t=-\mathcal{K}(\vartheta )\nabla \vartheta ,$$

$$\rho c_{\mathrm{v}}(\vartheta ){\vartheta }_t+\nabla \cdot \mathbf{q}=0,$$

$$\mathbf{q}+\tau (\vartheta ){\mathbf{q}}_t=-\mathcal{K}(\vartheta )\nabla \vartheta ;$$

$${\tau }_{\mathrm{GFS}}(\vartheta ):=\alpha \mathcal{K}(\vartheta )\vartheta ,$$

$$(\alpha \vartheta {\mathbf{q}}_t+\nabla \vartheta )\cdot \mathbf{q}\le 0,$$

$$\tau (\vartheta )=\frac{\mathcal{K}(\vartheta )(A_0+B_0{\vartheta }^n)}{\rho c_{\mathrm{v}}(\vartheta )},$$

$$q+\alpha K\vartheta q_t=-K{\vartheta }_x,$$

$$\rho c_{\mathrm{p}}{\vartheta }_t+q_x=0.$$

$${\left(\begin{array}{c}\vartheta \\ q\end{array}\right)}_t+\mathsf{A}{\left(\begin{array}{c}\vartheta \\ q\end{array}\right)}_x=-(\alpha K\vartheta {)}^{-1}\left(\begin{array}{c}0\\ q\end{array}\right),\text{where}\mathsf{A}=\left(\begin{array}{cc}0& \kappa /K\\ 1/(\alpha \vartheta )& 0\end{array}\right).$$

$$U(\vartheta )=\sqrt{\frac{\kappa }{\alpha K\vartheta }},$$

$$U_{\mathrm{wnl}}(\vartheta )=\sqrt{\frac{\kappa }{\alpha K{\vartheta }_0}}\left[1-\frac{1}{2}\left(\frac{\vartheta -{\vartheta }_0}{{\vartheta }_0}\right)\right](|\vartheta -{\vartheta }_0|\ll {\vartheta }_0),$$

$$U_0=\sqrt{\frac{\kappa }{\alpha K{\vartheta }_0}},$$

$$\theta =\vartheta /{\vartheta }_0,x^{\circ }=x/L,t^{\circ }=t(\kappa /L^2),q^{\circ }=q(L/(K{\vartheta }_0)),$$

$${\left(\begin{array}{c}\theta \\ q\end{array}\right)}_t+\left(\begin{array}{cc}0& 1\\ 1/(\lambda \theta )& 0\end{array}\right){\left(\begin{array}{c}\theta \\ q\end{array}\right)}_x=-(\lambda \theta {)}^{-1}\left(\begin{array}{c}0\\ q\end{array}\right).$$

$$\lambda :=\alpha \kappa K{\vartheta }_0/L^2,$$

$$\theta (0,t)=1+H_{\mathrm{p}}(t,t_{\mathrm{w}})\psi (t),\text{where}{H}_{\mathrm{p}}(t,t_{\mathrm{w}}):=H(t)-H(t-t_{\mathrm{w}}),$$

$$2\frac{\mathfrak{d}a}{\mathfrak{d}t}=-a({\lambda }^{-1}+a).$$

$$\frac{\mathfrak{d}[[\mathfrak{F}]]}{\mathfrak{d}t}=[[{\mathfrak{F}}_t]]+c_0[[{\mathfrak{F}}_z]],$$

$$a(t)=-|{\alpha }^{\ast }|{\left\{1-\left[1+\frac{|{\alpha }^{\ast }|}{a(0)}\right]\exp (\frac{1}{2}t/\lambda )\right\}}^{-1},$$

$${\alpha }^{\ast }:=-{\lambda }^{-1}.$$

$$t_{\infty }=2\lambda \ln \left[\frac{a(0)}{a(0)+|{\alpha }^{\ast }|}\right](0<t_{\infty }<\infty ).$$

$$V_t+c_0^{-2}(1+V_x)V_{tt}-V_{xx}=0,(x,t)\in (0,1)\times (-\infty ,t_{\mathrm{r}});$$

$$V_x(0,t)=\delta H_{\mathrm{p}}(t,t_{\mathrm{r}})\sin (\pi t),V_x(1,t)=0,t\in (-\infty ,t_{\mathrm{r}});$$

$$V(x,0)=0,V_t(x,0)=0,x\in (0,1).$$

$$V:\Omega \subset {\mathbb{R}}^2\mapsto \mathbb{R}s.t.V_x=-1+\theta ,$$

$$[[V_{xx}]]=-c_0^{-1}[[V_{xt}]]=\frac{|{\alpha }^{\ast }|}{c_0}{\left\{1-\left[1+\frac{|{\alpha }^{\ast }|}{\delta \pi }\right]\exp (\frac{1}{2}{c}_0^{2}t)\right\}}^{-1},$$

$$t_{\infty }=2\lambda \ln \left[\frac{\delta \pi }{\delta \pi +|{\alpha }^{\ast }|}\right].$$

$$\begin{gathered} \frac{V_m^{k+1}-V_m^{k-1}}{2(\Delta t)}+c_0^{-2}\left(1+\frac{V_{m+1}^k-V_{m-1}^k}{2(\Delta x)}\right)\left(\frac{V_m^{k+1}-2V_m^k+V_m^{k-1}}{(\Delta t{)}^2}\right) \\ -\frac{V_{m+1}^k-2V_m^k+V_{m-1}^k}{(\Delta x{)}^2}=0, \end{gathered}$$

$$V_{m}^{k+1}=\Bigg{\{}\tfrac{1}{2}(\Delta t)+c_{0}^{-2}\Bigg{[}1+\frac{V_{m+1}^{k}-V_{m-1}^{k}}{2(\Delta x)}\Bigg{]}\Bigg{\}}^{-1}\\ \times\Bigg{\{}R^{2}(V_{m+1}^{k}-2V_{m}^{k}+V_{m-1}^{k})+\tfrac{1}{2}(\Delta t)V_{m}^{k-1}\\ +c_{0}^{-2}(2V_{m}^{k}-V_{m}^{k-1})\Bigg{[}1+\frac{V_{m+1}^{k}-V_{m-1}^{k}}{2(\Delta x)}\Bigg{]}\Bigg{\}},$$

$$V_{-1}^k=V_1^k-2\delta (\Delta x)\sin (\pi t_k),V_{M+1}^k=V_{M-1}^k(k=1,2,3,\dots ,N),$$

$$V_m^0=0,V_m^1=V_m^0(m=-1,0,1,\dots ,M+1).$$

$$h-v\alpha Kf{h}^{\prime }=-K{f}^{\prime },$$

Full text

keyword Graffi–Franchi–Straughan conductors; Maxwell–Cattaneo law; temperature-rate waves; traveling wave solutions; Lambert $W$-function

1 Introduction

Absent the presence of source terms, the most general form of the system that describes the propagation of second-sound (i.e., thermal waves) in the ‘class’ of rigid, homogeneous and isotropic, solids that the present authors have termed Graffi–Franchi–Straughan (GFS) conductors [5] reads

[TABLE]

where $\vartheta(>0)$ denotes the absolute temperature, ${\bf q}$ is the heat flux vector, $c_{\rm v}(\vartheta)(>0)$ is the constant-volume specific heat [3], $\mathcal{K}(\vartheta)(>0)$ is the thermal conductivity, $\alpha$ is a positive constant that carries (SI) units of ${\rm m}\cdot{\rm sec}/{\rm W}$, and $\rho(>0)$ is the (constant) mass density of the conductor. Eq. (1.1a), we observe, is a special case of the Maxwell–Cattaneo (MC) law, i.e., a special case of the flux constitutive relation [11, 18, 19, 23, 24, 28]

[TABLE]

specifically, the former follows from the latter on setting $\tau(\vartheta)=\tau_{\rm GFS}(\vartheta)$, where

[TABLE]

which is the distinguishing feature of GFS conductors. Here, $\tau(\vartheta)(>0)$ denotes the relaxation time for phonon processes that are dissipative because they do not conserve phonon momentum [11].

Recalling arguments from an earlier (unpublished) contribution by Graffi, Franchi and Straughan [16], in 1994, used the derivation of Eq. (1.1a) to illustrate a theoretically-motivated means by which the temperature dependence of $\tau$ may arise. In particular, they noted that Eq. (1.1a) identically satisfies the following generalization of the Clausius–Duhem inequality:

[TABLE]

which Franchi and Straughan [16, p. 728] attribute to Graffi; see also Franchi [15] and Straughan [28, §1.2]. As alluded to above, however, GFS conductors must at present be regarded as hypothetical constructs since, to the best of our knowledge, the literature does not contain any examples of actual solids wherein q is described by Eq. (1.1a).

Nevertheless, GFS conductors exhibit a number of interesting mathematical properties that, from our perspective, make them worthy of investigation. For example, the GFS form of $\tau$ plays a critical role in establishing the following:

• •

The empirically based relation for second-sound in rigid solids

[TABLE]

where $n$, $A_{0}$, and $B_{0}$ are fitting parameters, has been shown to be applicable to NaF, for $10.0\leq\vartheta\leq 18.5\,$K, and Bi, for $1.4\leq\vartheta\leq 4.0\,$K; see, e.g., Refs. [10, 11] and those cited therein. On comparing with Sys. (1.1), it is easily seen that setting $n=1$, $A_{0}=0$, and $B_{0}=\alpha\rho c_{\rm v}(\vartheta)$ reduces Eq. (1.5) to Eq. (1.3). However, the fact that $B_{0}$ is a constant necessitates the additional requirement $c_{\rm v}(\vartheta):=$ const. This is true, as has long been known, in the case of many real solids when $\vartheta\gg\vartheta_{\rm D}$, where $\vartheta_{\rm D}$ denotes the Debye temperature of the solid in question; see, e.g., Ref. [3, §2]. From a strictly theoretical standpoint, then, Eq. (1.5) also applies to GFS conductors in the high-temperature regime, i.e., under conditions yielding $c_{\rm v}(\vartheta):=$ const.

• •

When $\tau(\vartheta)$ is given by Eq. (1.3) and $c_{\rm v}(\vartheta)$ and $\mathcal{K}(\vartheta)$ are both taken to be constant, the flux relation under Morro–Ruggeri (MR) theory [22] reduces to its simplest possible special case that still exhibits explicit dependence on $\vartheta_{t}$; i.e., the simplest possible special case in which Ref. [22, Eq. (48)] does not degenerate into a particular case of Eq. (1.2) (above).

• •

When $\tau(\vartheta)$ is given by Eq. (1.3), $\tilde{e}=\tilde{e}(\vartheta,{\bf q})$, the generalized expression for the internal energy density under Coleman–Fabrizio–Own (CFO) theory [9] (see also Refs. [10, 11, 22]), reduces to it simplest possible special case that still exhibits explicit dependence on ${\bf q}$; i.e., the simplest possible special case in which $\tilde{e}(\vartheta,{\bf q})$ does not degenerate into the classical expression for the internal energy density of a rigid solid.

The primary aim of this communication is to present numerical simulations of the second-sound phenomena that the present authors examined in Ref. [5] using only analytical methods. In particular, we simulate both temperature-rate waves and traveling waves predicted by the following special case of Sys. (1.1):

[TABLE]

As in Ref. [5], we have assumed $\mathcal{K}(\vartheta):=K$, where the constant $K$ is the value of the thermal conductivity at some reference temperature; we have taken $c_{\rm v}(\vartheta)$ to be a constant ($\Rightarrow\vartheta\gg\vartheta_{\rm D}$), but have also made use of the fact that $c_{\rm v}=c_{\rm p}$ under the rigid solid idealization222For details on the justification/rational behind the use of the approximation $c_{\rm v}\approx c_{\rm p}$ when modeling the flow of heat in real solids, see, e.g., Refs. [6, 14, 25]., where $c_{\rm p}$ denotes the constant-pressure specific heat; and we have confined our attention to one-dimensional (1D) heat flow along the $x$-axis, a propagation geometry which renders $\vartheta=\vartheta(x,t)$ and ${\bf q}=(q(x,t),0,0)$.

To this end, the present article is organized as follows. In Sect. 2, a review of the temperature-rate wave analysis carried out in Ref. [5] is presented. In Sect. 3, numerical simulations of temperature-rate waves are performed and results obtained are compared with our analytical findings. Then, in Sect. 4, a traveling wave analysis of Sys. (1.6) is performed and a two of the resulting solution profiles are studied numerically. And lastly, in Sect. 5, connections to other works are discussed and possible follow-on studies are noted.

2 Temperature-rate waves: Analytical results

2.1 Brief history and related works

By a temperature-rate wave333Also known as a temperature-rate discontinuity wave and a discontinuity (or acceleration) wave; see, e.g., Refs. [21] and [23, §7.4], respectively. we mean a singular surface, i.e., a wavefront, across which the first derivatives of the temperature field suffer a jump discontinuity; see, e.g., Refs. [21, 28]. What makes these waves so interesting is the fact that, under certain conditions, the jump amplitude can exhibit finite-time blow-up, even when the imposed thermal disturbance is continuous. Today, it is generally accepted that temperature-rate wave amplitude blow-up signals the formation of a thermal shock [28], i.e., a propagating surface across which $\vartheta$ itself suffers a jump.

2.2 Mathematical preliminaries: Characteristic speed

Letting $\kappa=K/(\rho c_{\rm p})$ denote the thermal diffusivity [14], we begin this subsection by recasting Sys. (1.6) in (equivalent) matrix form, specifically, as

[TABLE]

The eigenvalues, $\mu_{1,2}$, of the coefficient matrix satisfy the characteristic equation $\det(\mathsf{A}-\mu\mathsf{I}_{2})=0$, where $\mathsf{I}_{2}$ denotes the $2\times 2$ identity matrix; thus, $\mu_{1,2}=\pm U(\vartheta)$, where the characteristic speed of second-sound under Sys. (2.1) is

[TABLE]

and the characteristics of Sys. (2.1) are defined by ${\rm d}x/{\rm d}t=\pm U(\vartheta)$. Since $\mu_{1,2}\in\mathbb{R}$ and unequal, it follows that this (quasilinear) system is a strictly hyperbolic [20] one. Therefore, solutions of Sys. (2.1) satisfy the requirements of causality provided $U(\vartheta)$ is bounded, where from Eq. (2.2) we see that $U(\vartheta)\to\infty$ as $\vartheta\to 0$.

Since $\vartheta\gg\vartheta_{\rm D}$ is one of the assumptions on which Sys. (1.6) is based, the breakdown of our model as $\vartheta\to 0$ does not pose a difficulty for us. To gain deeper insight into the behavior of $U(\vartheta)$, it is instructive to consider finite- and small-amplitude thermal disturbances. Under the weakly-nonlinear and linear approximations Eq. (2.2) becomes

[TABLE]

where in this study $\vartheta_{0}$ denotes the initial temperature of the conductor, and

[TABLE]

respectively. Here, we see that $U_{\rm wnl}(\vartheta)>U_{0}$, when $\vartheta_{0}>\vartheta$, while $U_{\rm wnl}(\vartheta)<U_{0}$, when $\vartheta_{0}<\vartheta$. This, we observe, is the opposite of the behavior exhibited under the (bi-directional) model equations of classical acoustics; see, e.g., Ref. [17, §4($a$)], and note that $\wp$, the thermodynamic pressure in Ref. [17], corresponds to $\vartheta$.

Lastly, to simplify the forthcoming temperature-rate wave analysis, we now introduce the following non-dimensional variables:

[TABLE]

where $L$ denotes the conductor’s thickness, and recast Sys. (2.1) in non-dimensional form, viz.:

[TABLE]

where, for convenience, we have set

[TABLE]

all superscript circles have been omitted but should remain understood, and we note for later reference that $\mathcal{C}(\theta)=(\lambda\theta)^{-1/2}$ is the non-dimensional form of $U(\vartheta)$.

2.3 Formulation

Now consider a rigid conducting slab, whose (normalized) thickness is unity, wherein the temperature and heat flux are described by Sys. (2.6). We suppose the slab is stationary and that, initially, $q=0$ and the slab is at a uniform temperature $\theta(x,0)=1$ (i.e., $\vartheta(x,0)=\vartheta_{0}$) throughout. Beginning at time $t=0+$, let a temperature pulse of the form

[TABLE]

be applied to the boundary $x=0$, while the boundary $x=1$ is held at temperature $\theta(1,t)=1$. Here, the pulse duration (or width) $t_{\rm w}(>0)$ is a constant; the amplitude function $|\psi(t)|\in(0,1)$ is assumed to be continuously differentiable, nonzero on the interval $t\in(0,t_{\rm w})$, and such that $\psi(0)=0$ but $\psi_{t}(0)\neq 0$; and $H(\cdot)$ denotes the Heaviside unit step function.

Since Sys. (2.6) is strictly hyperbolic, a planar temperature wavefront $x=\Sigma(t)$, across which $[\![\theta]\!]=[\![q]\!]=0$, but $[\![\theta_{t}]\!]\neq 0$, begins propagating from the boundary $x=0$, along the positive $x$-axis, with speed $\mathcal{C}(\theta^{+})=\lambda^{-1/2}$ relative to the slab, where we observe that $\theta^{+}=1$ under the present formulation. Here, employing the standard notation of singular surface theory, $[\![\mathfrak{F}]\!]:=\mathfrak{F}^{-}-\mathfrak{F}^{+}$ denotes the amplitude of the jump in the value of the function $\mathfrak{F}=\mathfrak{F}(x,t)$ across $\Sigma(t)$, where $\mathfrak{F}^{\pm}\equiv\lim_{x\to\Sigma(t)^{\pm}}\mathfrak{F}(x,t)$ are assumed to exist, and $\pm$ superscripts correspond to the regions ahead of and behind $\Sigma$, respectively. Since, under this formulation, $\theta_{t}$ suffers a jump discontinuity across it, the surface $\Sigma$ is clearly a temperature-rate wave.

2.4 Amplitude evolution

Observing now that $[\![\theta_{t}]\!]$ is, at most, a function of only $t$, and referring the reader to Refs. [4, 28] for details, it is a straightforward matter to show that $\Sigma(t)=c_{0}t+x_{0}$, where we have set $c_{0}:=\lambda^{-1/2}$ for convenience and $x=x_{0}$ is the location of $\Sigma$ at $t=0$, and that the jump in $\theta_{t}$ satisfies the Bernoulli equation

[TABLE]

Here, use has been made of

[TABLE]

which is usually referred to as the kinematic condition of compatibility 444See Ref. [28, §4.1] and those cited therein; see also Bland [4, §6.9], who refers to this relation as ‘Hadamard’s lemma’., where $\mathfrak{d}/\mathfrak{d}t$, the 1D displacement derivative, gives the time-rate-of-change measured by an observer traveling with $\Sigma$; we have set $a(t):=[\![\theta_{t}]\!]$ for convenience; and it should be noted that, since the slab’s initial temperature was assumed to be constant, we took $\theta_{t}^{+}=0$ in deriving Eq. (2.9).

Making use of the substitution $a=1/\mathfrak{a}$, Eq. (2.9) is transformed into a linear ODE, which is easily integrated; its exact solution can be expressed as

[TABLE]

where the (negative) constant $\alpha^{*}$, known as the critical amplitude, is given by

[TABLE]

According to Eq. (2.11), $a(t)$ can evolve in any one of the following four ways:

1. (i)

If $a(0)>0$, then $a(t)\in(0,a(0))$ for $t>0$ and $a(t)\to 0$ from above as $t\to\infty$. 2. (ii)

If $a(0)<0$ and $|a(0)|<|\alpha^{*}|$, then $a(t)\in(a(0),0)$ for $t>0$ and $a(t)\to 0$ from below as $t\to\infty$. 3. (iii)

If $a(0)=\alpha^{*}$, then $a(t)=\alpha^{*}$ for all $t\geq 0$. 4. (iv)

If $a(0)<0$ and $|a(0)|>|\alpha^{*}|$, then $a(t)<a(0)$ for $t>0$ and, moreover, $|a(t)|\to\infty$ as $t\to t_{\infty}$, where

[TABLE]

2.5 Stability results

While we have obtained the exact solution of Eq. (2.9), it is nevertheless instructive to investigate the steady-state behavior of $a(t)$ using qualitative methods; i.e., to examine the stability characteristics of the equilibrium solutions $\bar{a}=\{0,-|\alpha^{*}|\}$, which of course correspond to the roots of the quadratic equation $-a(|\alpha^{*}|+a)=0$.

As a phase plane analysis reveals, $\bar{a}=-|\alpha^{*}|$ is always unstable while $\bar{a}=0$ is always stable. This means that a bifurcation does not occur in the case of Eq. (2.9); i.e., there is no interchange of stability between the two equilibria of this ODE. The instability of $\bar{a}=-|\alpha^{*}|$ also means that the constant solution in Case (iii) is unstable as well; i.e., any discrepancy, however small, in achieving $a(0)=\alpha^{*}$ will yield either Case (ii) or (iv).

3 Temperature-rate waves: Numerical results

While interesting and useful, temperature-rate wave results do not provide any information on the behavior of the temperature field behind $\Sigma$

Hence, to explore this aspect of the GFS model, and to illustrate the most important findings of Sect. 2.4, we now turn to computational methods. In this section, we present a series of numerical simulations based on the slab initial-boundary value problem (IBVP) formulated in Sect. 2.3, which we now express as:

[TABLE]

Here, so that Sys. (2.6) could be recast as a single PDE, we have introduced

[TABLE]

where $\Omega=\{(x,t):0<x<1,-\infty<t<t_{\rm r}\}$. Furthermore, $\psi(t)=\delta\sin(\pi t)$, where $|\delta|\in(0,1)$ is a constant; $t_{\rm w}=t_{\rm r}$, where $t_{\rm r}=1/c_{0}$ is the time required for $\Sigma$ to complete its initial transit of the interval $0<x<1$ (i.e., $\Sigma(t_{\rm r})=1$); and of course $x_{0}=0$.

In the case of IBVP (3.1) the temperature-rate wave amplitude expression, i.e., Eq. (2.11), becomes

[TABLE]

where the jump $[\![V_{xx}]\!]$ was determined using the expression for $[\![V_{xt}]\!]$, the fact that $[\![V_{x}]\!]=0$, and the $[\![\mathfrak{F}]\!]=0$ special case of Eq. (2.10), while the expression for the blow-up time, i.e., Eq. (2.13), assumes the form

[TABLE]

Also, it should be noted that the amplitude of the jump in the time derivative of the boundary condition at $x=0$ across the plane $t=0$ is $[\![V_{xt}]\!]|_{t=0}=a(0)=\delta\pi$.

Since an exact analytical solution does not appear to be possible, we turn to the calculus of finite differences and introduce the mesh points $(x_{m},t_{k})$, where $x_{m}=m(\Delta x)$ for each $m=-1,0,1,\ldots,M+1$ and $t_{k}=k(\Delta t)$ for each $k=0,1,2,\ldots,N$. Here, the spatial- and temporal-step sizes are defined as $\Delta x=1/M$ and $\Delta t=T/N$, respectively, where $M(\geq 2)$ and $N(\geq 2)$ are integers and $T\in(0,t_{\rm r})$ is the right-hand endpoint555Of course, if $t_{\infty}\in(0,t_{\rm r})$, then the restriction on $T$ becomes $T\in(0,t_{\infty})$. of the temporal interval over which the solution of IBVP (3.1) shall be computed.

With our (1D) mesh established, and guided by the treatment of similar equations presented in Refs. [2, 17], we construct the following simple discretization of Eq. (3.1a):

[TABLE]

where $V_{m}^{k}\approx V(x_{m},t_{k})$. On setting $R=(\Delta t)/(\Delta x)$ and then solving for $V_{m}^{k+1}$, the most advanced time-step approximation, we obtain the (explicit) finite difference scheme (FDS)

[TABLE]

which holds for each $m=0,1,2,\ldots,M$ and $k=1,2,3\ldots,N-1$. In turn, discretization of the boundary conditions gives

[TABLE]

where we note our use of the ghost points666A numerical device that allows us to discretize the Neumann boundary conditions of our problem using centered-difference quotations, i.e., consistent with how the spatial derivatives in Eq. (3.1a) are discretized in our finite difference scheme; see, e.g., Ref. [29]. $m=-1,M+1$, while the initial conditions become

[TABLE]

In Figs. 1–3 we have presented temperature profile plots corresponding to Cases (i), (iii), and (iv), respectively. These time-sequence plots depict the evolution of the $V_{x}$ vs. $x$ solution profile under IBVP (3.1), with $V_{x}$ normalized by $\delta$, during $\Sigma$’s initial transit of the slab. The curves shown in solid black were produced from data sets computed by a simple algorithm which implemented FDS (3.6) on a desktop computer running Mathematica (ver. 11.2). Interpolations between the points were then accomplished using the cubic interpolation routine that is a built-in part of this software package. The red broken lines, which were generated from Eq. (3.3), have been included to illustrate the behavior of the temperature-rate wave amplitudes; as $V_{xx}^{+}=0$ under IBVP (3.1), the slopes of these lines give the values of $V_{xx}^{-}$, at their points of tangency to the solution profiles, at the indicated times.

For consistency across these three figures, and ease of computation, we have selected the common values $\lambda=1.2$ ($\Rightarrow c_{0}\approx 0.9129$, $t_{\rm r}\approx 1.0954$) and $M=2500$, $N=5000$, $T=1$ ($\Rightarrow R=1/2$). It should be noted that of those we tested, with $R=1/2$ fixed, $\lambda=1.2$ was the smallest value of $\lambda$ for which both FDS (3.6) was numerically stable and we could place, in the case of Fig. 3, $x_{\infty}=c_{0}t_{\infty}$ very close to, but to the left of, the boundary $x=1$. The values of $M$ and $N$ selected were based on a heuristic search to find the smallest such values, subject to $R=1/2$, that accurately captured the manifestation of the temperature-rate wave on the temperature profile in the last frame of Fig. 3.

In Fig. 1 we observe, as predicted in Case (i), the slope of the profile at the wavefront decreasing to zero, as $t\to\infty$, when $\delta>0$. In Fig. 2 we see, as predicted in Case (iii), the slope of the profile at the wavefront remaining constant; specifically, $V_{xx}^{-}=[\![V_{xx}]\!]=-1/\lambda$. In contrast, Fig. 3777. Note that in Fig. 3, $T<t_{\infty}<t_{\rm r}$., which captures approximately $87\%$ of the ‘lifetime’ of $\Sigma$, clearly illustrates the exponential increase in $|[\![V_{xx}]\!]|$ as $t\to t_{\infty}$, as predicated in Case (iv). In particular, the last frame of Fig. 3 shows the slope on the leading side of our solution profile becoming nearly vertical at the wavefront, strongly suggesting that a thermal shock is about to form.

4 Traveling wave analysis888The reader

should be aware that the analysis presented in Ref. [5, §3] contains a number of omissions, misstatements, and misprints. In the present section, these issues have, without identification nor comment, all been remedied/corrected.

4.1 Associated ODE, jump magnitude

Assuming right-running waveforms propagating along the $x$-axis, we take the dependence of $\vartheta$ and $q$ on $x$ and $t$ to be of the form $\vartheta(x,t)=f(\xi)$ and $q(x,t)=h(\xi)$, where $\xi:=x-vt$ is the wave variable and $v(>0)$ is the (constant) wave speed. On substituting these ansatzs into Sys. (1.6) we obtain, after simplifying, the system of ODEs

[TABLE]

a (trivial) solution of which, we observe, is

[TABLE]

Here, a prime denotes $\textrm{d}/\textrm{d}\xi$ and $\vartheta^{\bullet}(>0)$ is a constant.

Now eliminating $h$ between the equations of Sys. (4.1), after integrating Eq. (4.1b) once, and then assuming999That is, we are seeking kink [1], and kink-like, traveling wave solutions. that $f(\xi)\to\vartheta_{\rm r}$, $f^{\prime}(\xi)\to 0$ as $\xi\to-\infty$, we obtain the following Abel equation [13] for the temperature field:

[TABLE]

which is the associated ODE of Sys. (1.6). Here, we recall that $\kappa=K/(\rho c_{\rm p})$ is the thermal diffusivity; the constant $\vartheta_{\rm r}$ denotes a reference state value of $\vartheta$; and enforcement of the asymptotic condition gives $\mathfrak{K}_{1}=-\rho c_{\rm p}v\vartheta_{\rm r}$, where $\mathfrak{K}_{1}$ is the resulting constant of integration.

An inspection of Eq. (4.3) reveals that $\bar{f}=\vartheta_{\rm r}$, the only equilibrium point of this ODE, is unstable for $v>v_{\rm a}$, but stable for $v<v_{\rm a}$, where

[TABLE]

Here, we observe that $v=v_{\rm a}$ is a degenerate case in the following sense: $f=\vartheta_{\rm r}$ satisfies Eq. (4.3), and it is also true that

[TABLE]

i.e., taking $v=v_{\rm a}$ causes $f^{\prime}$ to exhibit a jump discontinuity at $\bar{f}=\vartheta_{\rm r}$, the magnitude of which is

[TABLE]

Introducing now the dimensionless temperature $\mathcal{T}=f/\vartheta_{\rm r}$ and dimensionless wave variable $\eta=\xi/\ell$, Eq. (4.3) is reduced to

[TABLE]

Here, $\ell(>0)$ is a characteristic length; we have set

[TABLE]

where $\sigma=1$ implies $v=v_{\rm a}$; and $c$, the dimensionless version of the wave speed $v$, is given by

[TABLE]

where we note that $c_{\rm a}=v_{\rm a}\ell/\kappa$ (i.e., $c_{\rm a}$ is the dimensionless version of $v_{\rm a}$). Also, for later reference we observe that, in terms of the present dimensionless quantities, Eqs. (4.2) and (4.6) become

[TABLE]

where $\mathfrak{h}$ denotes the dimensionless version of $h$ and $\mathcal{T}^{\bullet}=\vartheta^{\bullet}/\vartheta_{\rm r}$, and

[TABLE]

respectively.

4.2 Complete stability results

A full phase plane analysis of Eq. (4.7) reveals the following:

1. (I)

If $\sigma>1$ and $\mathcal{T}_{\rm w}>1$, then $\overline{\mathcal{T}}=1$ is unstable (above) and $1<\mathcal{T}(\eta)<\infty$, where $\mathcal{T}(\eta)\to\infty$ as $\eta\to\infty$. 2. (II)

If $\sigma>1$ and $\sigma^{-1}<\mathcal{T}_{\rm w}<1$, then $\overline{\mathcal{T}}=1$ is unstable (below) and $\sigma^{-1}\leq\mathcal{T}(\eta)<1$. 3. (III)

If $\sigma>1$ and $0<\mathcal{T}_{\rm w}<\sigma^{-1}$, then $-\infty<\mathcal{T}(\eta)\leq\sigma^{-1}$, where $\mathcal{T}(\eta)\to-\infty$ as $\eta\to-\infty$. 4. (IV)

If $\sigma=1$ and $\mathcal{T}_{\rm w}>1$, then $1\leq\mathcal{T}(\eta)<\infty$, with $\mathcal{T}(\eta)\to\infty$ as $\eta\to\infty$. 5. (V)

If $\sigma<1$ and $\mathcal{T}_{\rm w}>\sigma^{-1}$, then $\sigma^{-1}\leq\mathcal{T}(\eta)<+\infty$, where $\mathcal{T}(\eta)\to+\infty$ as $\eta\to+\infty$. 6. (VI)

If $\sigma<1$ and $1<\mathcal{T}_{\rm w}<\sigma^{-1}$, then $\overline{\mathcal{T}}=1$ is stable (above) and $1<\mathcal{T}(\eta)\leq\sigma^{-1}$. 7. (VII)

If $\sigma<1$ and $0<\mathcal{T}_{\rm w}<1$, then $\overline{\mathcal{T}}=1$ is stable (below) and $-\infty<\mathcal{T}(\eta)<1$, with $\mathcal{T}(\eta)\to-\infty,1$ as $\eta\to\mp\infty$, respectively.

Here, $\mathcal{T}(0)=\mathcal{T}_{\rm w}$, where $\mathcal{T}_{\rm w}$ is a (known) positive constant; however, as Cases (I), (IV), and (VI) shall be of particular interest to us, we hereafter limit our attention to $\mathcal{T}_{\rm w}>1$.

Returning to Eq. (4.7), we separate variables and integrate; this yields, after then applying and enforcing the condition at $\eta=0$,

[TABLE]

where $\mathfrak{K}_{2}(\mathcal{T}_{\rm w})=\sigma(\mathcal{T}_{\rm w}-1)+(\sigma-1)\ln(\mathcal{T}_{\rm w}-1)$.

In the next two subsections, the cases of $\sigma\geq 1$ and $\sigma\in(0,1)$ are treated consecutively.

4.3 The case $\sigma\geq 1$

For $\sigma\geq 1$, the integral curves take the form

[TABLE]

where $W_{0}(\cdot)$ denotes the principal branch of the Lambert $W$- function [12] and we have set $\eta_{\rm c}:=c_{\rm a}^{-1}(1-\mathcal{T}_{\rm w})$. For the case $\sigma=1$, the solution profile is seen to be a piecewise-linear, but continuous, function of $\eta$; it was constructed by joining together, at the point $\eta=\eta_{\rm c}$, Eq. (4.10), with $\mathcal{T}^{\bullet}=1$, and the $\sigma=1$ special case of Eq. (4.12).

From Eq. (4.13) we find that $\mathcal{T}\!\to\!\infty$ as $\eta\!\to\!\infty$. It can, however, be shown that the temperature gradient in this case is bounded; specifically, $0\leq\textrm{d}\mathcal{T}(\eta)/\textrm{d}\eta\leq c_{\rm a}$, for all $\sigma\geq 1$, where in this subsection the temperature gradient is given by

[TABLE]

From Eq. (4.14) it is clear that $\sigma=1$ corresponds to Case (iii), i.e., the constant amplitude case, of our temperature-rate wave analysis (recall Sect. 2.4). The sequence shown in Fig. 4 depicts the steepening of the temperature gradient profile as $\sigma\to 1$ (from above); in this limit, the $\textrm{d}\mathcal{T}(\eta)/\textrm{d}\eta$ vs. $\eta$ profile tends to a step function, i.e., a temperature-rate wave in the present setting, of (jump) magnitude $c_{\rm a}$.

4.4 The case $\sigma\in(0,1)$

On joining, at the point $\eta=\eta^{*}$ (see below), the $\sigma\in(0,1)$ special case of Eq. (4.12) to the constant-temperature solution given in Eq. (4.10), the piecewise-defined integral curve corresponding to this range of $\sigma$-values is readily constructed, viz.:

[TABLE]

Here, we note the restriction $\mathcal{T}_{\rm w}\in(1,\sigma^{-1})$ and observe that the value of $\mathcal{T}(\eta)$ on the interval $\eta\leq\eta^{*}$ follows on setting $\mathcal{T}^{\bullet}=1/\sigma$ in Eq. (4.10). Also, $\eta^{*}(<0)$ is given by

[TABLE]

where Eq. (4.16) was obtained by setting the argument of $W_{0}$ (in Eq. (4.15)) equal to $-1/\textrm{e}$ and then solving for $\eta$, where $-1/\textrm{e}$ is a branch point of the $W$-function; again, see Ref. [12].

Along with the fact that Eq. (4.15) describes a bounded, continuous, waveform, for which $\overline{\mathcal{T}}=1$ is a stable equilibrium, Fig. 5 also reveals that for $0<\sigma<1$, the GFS traveling wave profile exhibits what Zel’dovich and Raizer [31, Fig. 10.3b] refer to as a (preheating) ‘tongue’.

5 Closure

From the mathematical standpoint, the present analysis has shown that the qualitative behavior of Sys. (1.6) is very much like that exhibited by the 1D version of the ‘Darcy–Jordan’ model101010Also known as the ‘Jordan–Darcy’ model; see, e.g., Ref. [26] and those cited therein. [17] (of poroacoustics) when the fluid phase consists of a retrograde fluid111111Fluids that exhibit rarefaction (or ‘negative’) shocks; see, e.g., Ref. [30] and those cited therein. Note also that retrograde fluids correspond to $\beta<0$ in Ref. [17], wherein $\beta$ denotes the coefficient of nonlinearity.. This analogy is most evident in the case of temperature-rate waves; recall the behavior of $U_{\rm wnl}(\vartheta)$ (see Sect. 2.2), and observe that $-V_{x}$ corresponds to $p^{\prime}$, where $p^{\prime}$ is used in Ref. [17, §4] to denote the dimensionless over pressure in a regular fluid.

It is also noteworthy that the behavior of the waveforms observed in Sect. 4 under the $\sigma\geq 1$ and $\sigma\in(0,1)$ cases corresponds to taking $\epsilon>0$ and $\epsilon<0$, respectively, in Christov and Jordan [8, pp. 1126–1128], who examined traveling waves under the MC law with $\mathcal{K}(\vartheta)$ a linear function of $\vartheta$ and $\tau(\vartheta):=$ const.

With regard to possible follow-on studies, the obvious next step from the numerical standpoint is to examine thermal shock phenomena under Sys. (1.6) using what are known as ‘shock capturing’ schemes (see, e.g., Ref. [8] and those cited therein), which are more elaborate than the simple explicit scheme we employed in Sect. 3. From the analytical standpoint, extensions of the present study might include performing the above analyzes on the GFS special cases of MR theory and CFO theory; recall the second and third bulleted items in Sect. 1.

On the other hand, Sys. (1.1) could be recast in terms of what Straughan [27] has termed ‘Cattaneo–Christov’ theory. Under this generalization of the MC law, which Christov [7] proposed in 2009, the simple partial time derivative that acts on ${\bf q}$ would be replaced by a Lie derivative—one corresponding to Oldroyd’s upper convected derivative—which would then make Sys. (1.1) applicable to moving (GFS) conductors; to this latter point, see also the first footnote in Ref. [9, p. 136].

Acknowledgments

The authors are grateful to the anonymous referee for his/her helpful comments and for bringing Ref. [10] to their attention. S.C. thanks the financial support of G.N.F.M.–I.N.d.A.M., I.N.F.N. and Università di Roma La Sapienza, Rome, Italy.