*2.3. Dual-Phase-Lag Approximation Scheme for Test Structure*

In order to obtain the solution of temperature distribution problem inside the investigated cross-sectional area of the structure, the Dirichlet boundary conditions were considered at the bottom, while left, right and the top parts have been modeled using Neumann boundary conditions. Described situation can be characterized by the following equations:

$$\begin{aligned} T\_k(t) &= 0, \quad &t \in \mathbb{R}\_+ \cup \{0\}, \\ &k \in \{1, 2, 3, \dots, n\_x\} \end{aligned} \tag{12}$$

$$\begin{array}{c} \mathbf{q}\_{\mathbf{k}}(t) = \mathbf{0}, \ t \in \mathbb{R}\_{+} \cup \{0\}, \\ k \in \left\{ n\_{\mathbf{x}} + 1, 2 \cdot n\_{\mathbf{x}} + 1, \dots, \left( n\_{\mathbf{y}} - 1 \right) \cdot n\_{\mathbf{x}} + 1 \right\} \cup \left\{ n\_{\mathbf{x}}, 2 \cdot n\_{\mathbf{x}}, \dots, \left( n\_{\mathbf{y}} - 1 \right) \cdot n\_{\mathbf{x}} \right\} \cup \\ \cup \left[ \left( n\_{\mathbf{y}} - 1 \right) \cdot n\_{\mathbf{x}} + 1, \left( n\_{\mathbf{y}} - 1 \right) \cdot n\_{\mathbf{x}} + 2, \dots, n\_{\mathbf{y}} \cdot n\_{\mathbf{x}} \right] \end{array} \tag{13}$$

*Energies* **2020**, *13*, 2379

Considering imposed boundary conditions and the FDM assumptions, the DPL model can be explained by the following system:

$$\begin{cases} \mathbf{M} \cdot \mathbf{\dot{T}}(t) + \mathbf{D} \cdot \mathbf{\dot{T}}(t) + \mathbf{K} \cdot \mathbf{T}(t) = \mathbf{b} \cdot u(t) \\ \mathbf{y}(t) = \mathbf{c}^T \cdot \mathbf{T}(t) \end{cases} \quad t \in \mathbb{R}\_+ \cup \{0\} \tag{14}$$

where index *T* means a transposition, while **T**, **T** ˙ , and **T** ¨ are *nx*·*ny*-element vectors reflecting the temperature function and its first and the second time derivatives, respectively. Moreover, taking into account the way of nodes numbering, matrices **Idiag**, **MFDM**, **M**, **D**, **K**, and **c***T*; vectors **cv**, **k**, <sup>τ</sup>**q**, τ**T**, **b**, and **y**; and the function *u*(*t*) may be reflected as indicated below:

**Idiag**, **MFDM**,**M**, **D**, **K**, **c***T* ∈ *Rnx*·*ny* <sup>×</sup>*nx*·*ny* , **cv**, **k**, <sup>τ</sup>**q**, τ**T b**, **y** ∈ *Rnx*·*ny* × 1, *u* ∈ *R* **MFDM** = ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ −31 0 ... 100 ... 000 1 −4 1 ... 010 ... 000 0 1 −4 ... 001 ... 000 ... ... ... ... ... ... ... ... ... ... ... 100 ... −41 0 ... 100 010 ... 1 −4 1 ... 010 001 ... 0 1 −4 ... 001 ... ... ... ... ... ... ... ... ... ... ... 000 ... 100 ... −31 0 000 ... 010 ... 1 −3 1 000 ... 001 ... 0 1 −2 ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ (15) **M** = diag<sup>τ</sup>**q** ◦ diag(**cv**) ◦ **Idiag** (16) 1

$$\mathbf{D} = \text{diag}(\mathbf{c}\_{\mathbf{v}}) \circ \mathbf{I}\_{\text{diag}} - \frac{1}{\left(\Delta x\right)^{2}} \cdot \text{repmat}(\mathbf{r}\mathbf{r}) \circ \text{repmat}(\mathbf{k}) \circ \mathbf{M} \text{FDM} \tag{17}$$

$$\mathbf{K} = -\frac{1}{\left(\Delta x\right)^{2}} \mathbf{repmat}(\mathbf{k}) \circ \mathbf{M} \mathbf{FDM} \tag{18}$$

$$\mathbf{b} = a \cdot \begin{bmatrix} 0 & \cdots & 0 & 1 & \cdots & 1 & 0 & \cdots & 0 \end{bmatrix}^T, \qquad a \in \mathbb{R}\_+ \tag{19}$$

$$\mathbf{c} = \mathbf{I}\_{\text{diag}} \tag{20}$$

where **Idiag** is an identity matrix, operator ◦ reflects matrices multiplication resulting in a Hadamard product, and matrix function diag (·) creates a diagonal matrix from a vector, while Matlab repmat function replicates a given vector and composes a matrix of required dimensions. It is also worth highlighting that non-zero values in vector **b** are observed only in the case of mesh nodes characterizing a heating area.

˙

The system in Equation (14) of differential equations of the second order was equivalently transformed into the following system of the linear first-order equations [39]:

$$\begin{cases} \mathbf{E} \cdot \dot{\overline{\mathbf{T}}}(t) = \mathbf{A} \cdot \overline{\mathbf{T}}(t) + \mathbf{B} \cdot u(t) \\ \mathbf{y}(t) = \mathbf{c}^T \cdot \overline{\mathbf{T}}(t) \end{cases} \qquad t \in \mathbb{R}\_+ \cup \{0\} \tag{21}$$

where **T** and **T** are <sup>2</sup>·*nx*·*ny*-element vectors. The first part of vector **T** consists of elements of vector **T**, while its second part includes elements of **T** ˙ . On the other hand, the first and the second parts of the vector ˙ **T** consists of elements of vectors **T** ˙ and **T** ¨ , respectively. Moreover, matrices **E**, **A**, **B**, and **C***<sup>T</sup>* can be represented in the following way [39,40]:

$$\mathbf{E} = \begin{bmatrix} \mathbf{I}\_{\text{diag}} & \boldsymbol{\Theta} \\ \boldsymbol{\Theta} & \mathbf{M} \end{bmatrix} \cdot \mathbf{A} = \begin{bmatrix} \boldsymbol{\Theta} & \mathbf{I}\_{\text{diag}} \\ -\mathbf{K} & -\mathbf{D} \end{bmatrix} \cdot \mathbf{B} = \begin{bmatrix} \boldsymbol{\Theta}\_{1} \\ \mathbf{b} \end{bmatrix} \cdot \mathbf{C} = \begin{bmatrix} \mathbf{c} & \boldsymbol{\Theta} \\ \boldsymbol{\Theta} & \mathbf{c} \end{bmatrix} \tag{22}$$

where Θ and Θ1 are null matrices. Moreover, we ge<sup>t</sup> the following:

**E**, **A**, **C***<sup>T</sup>* ∈ *R*2·*nx*·*ny* <sup>×</sup>2·*nx*·*ny* , **B** ∈ *R*2·*nx*·*ny* ×1, **I**, Θ ∈ *Rnx*·*ny* <sup>×</sup>*nx*·*ny* , Θ1 ∈ *Rnx*·*ny* ×1

Thus, the full authors' approximation scheme for the DPL model in two-dimensional Euclidean space is described by the system in Equation (21), including explanations formulated in Equations (15)–(22) and boundary conditions in Equations (12) and (13). The final solution of temperature-distribution changes over time in the analyzed cross-section of the test structure, based on the prepared approximation scheme, is determined by using the BDF method of a variable order between 1 and 5.
