Analysis of Algorithm Efficiency for Heat Diffusion at Nanoscale Based on a MEMS Structure Investigation

Abstract

This paper presents an analysis of the time complexity of algorithms prepared for solving heat transfer problems at nanoscale. The first algorithm uses the classic Dual-Phase-Lag model, whereas the second algorithm employs a reduced version of the model obtained using a Krylov subspace method. This manuscript includes a description of the finite difference method approximation prepared for analysis of the real microelectromechanical system (MEMS) structure manufactured by the Polish Institute of Electron Technology. In addition, an approximation scheme of the model, as well as the Krylov subspace-based model order reduction technique are also described. The paper considers simulation results obtained using both investigated versions of the Dual-Phase-Lag model. Moreover, the relative error generated by the reduced model, as well as the computational complexity of both algorithms, and a convergence of the proposed approach are analyzed. Finally, all analyses are discussed in detail.

1. Classical and Modern Heat Transfer Description

Currently, deep development in technology has increased interest in heat transfer modeling. New thermal problems have been observed due to the continuous reduction of the size of electronic devices or miniaturization of integrated circuits. For instance, in such small electronic appliances, a significant rise of operational frequency and a rapid increase of generated internal heat density have been noticed which have a meaningful influence on temperature rise during the operation of the device.

It is worthwhile to highlight that it is extremely important to ensure the appropriate cooling conditions, however, this issue is very problematic in the case of nanosized electronic devices and can result in unstable and improper operation of the device. Moreover, most of the damages and malfunctions of the device are caused by unsuitable operation and thermal problems. Thus, thermal analysis is very important and is one of the most crucial steps in the design and planning of modern electronic appliances.

The heat transfer problems have been modeled using Fourier’s theory [1]. This method is based on Fourier’s law and the Fourier–Kirchhoff (FK) equation, and can be expressed as follows:

$$\{\begin{array}{cc}\begin{array}{c}\mathit{q}\left(\mathit{x},t\right)=-k\cdot \nabla T\left(\mathit{x},t\right)\\ c_v\frac{\partial T\left(\mathit{x},t\right)}{\partial t}+\nabla \circ \mathit{q}\left(\mathit{x},t\right)=q_V\left(\mathit{x},t\right)\end{array}& \mathit{x}\in R^n,n\in N,t\in R_{+}{{\displaystyle \cup }}^{}\left\{0\right\}\end{array}$$

(1)

In the above system of equations, the heat flux density vector is represented by q, thermal conductivity of analyzed material is expressed by k, and T describes the function of temperature rise. Moreover, the volumetric heat capacity and value of internally generated heat are represented by c_v and q_V, respectively.

The Fourier–Kirchhoff method has been successfully applied for two centuries. However, in modern nanosized structures, application of the FK method has some serious limitations [2,3,4]. First, the assumption of infinite speed propagation of the heat is postulated. Moreover, the instantaneous change of heat flux or temperature gradient should also be taken into account as a non-physical behavior of the structures. The phenomena mentioned have not been empirically proven in the case of nanosized electronic structures [5,6]. Thus, a new methodology, which significantly improves the FK method in the case of modern electronic structures, should be established.

In the mid-1990s, a new approach, called Dual-Phase-Lag (DPL), was introduced by Tzou as a more appropriate choice for modeling temperature changes in a nanosized electronic structure [7]. The DPL model based on the FK theory, however, as previously mentioned, includes improvements such as time lags [8]. These lags express the needed time change in the heat flux density, as well as the temperature gradient.

The lags mentioned above are represented by different values expressed by τ_q, which is related to heat flux time lag, and τ_T which represents the temperature time lag. Taking into consideration new time lags, the DPL model can be described by the following equations:

$$\begin{array}{cc}\{\begin{array}{c}\nabla \circ \mathit{q}\left(\mathit{x},t\right)=-c_v\frac{\partial T\left(\mathit{x},t\right)}{\partial t}+q_V\left(\mathit{x},t\right)\\ k{\tau }_T\frac{\partial \nabla T\left(\mathit{x},t\right)}{\partial t}+{\tau }_q\frac{\partial \mathit{q}\left(\mathit{x},t\right)}{\partial t}=-k\cdot \nabla T\left(\mathit{x},t\right)-\mathit{q}\left(\mathit{x},t\right)\end{array}& \mathit{x}\in R^n,n\in N,t\in R_{+}{{\displaystyle \cup }}^{}\left\{0\right\}\end{array}$$

(2)

The DPL model can be successfully applied instead of the classical FK model due to the fact that it is appropriate for parabolic partial differential equations, as well as for hyperbolic equations (see also [9,10]).

However, some disadvantages have also appeared. The DPL model is classified as a model with greater complexity than the FK model. Thus, to carry out the simulation, a longer computation time is needed, especially for complex electronic structures.

Considering the above limitation of the DPL model, the main aim of this paper is to implement a DPL-based method which reduces the time of simulation and can be as accurate as the original DPL model. One approach which can be applied is the Krylov subspace-based model order reduction method [11]. This approach significantly reduces the number of equations in the system describing an analyzed heat transfer problem.

The Krylov subspace method is used for the order reduction of large-scale systems, especially in a mechanical domain (e.g., J. White and J. G. Korvink papers, as well as [12,13]). However, the research presented in this paper is focused on aspects other than those of previously published papers. First of all, this manuscript includes the first application of the Krylov subspace-based method for the DPL heat transfer equation. Moreover, in previous research, the spectrum of mechanical systems generally contains frequencies from 0 to hundreds of Hz (sometimes up to 1–5 kHz). The model is validated typically for longer times, for example, 10 ms–1 s. However, our research includes significantly greater frequency values resulting in meaningfully shorter times and smaller DPL time lags, even hundreds of femtoseconds or a few nanoseconds. Therefore, the range of applicability of the described model order reduction methodology, presented in the manuscript, is different than other previous applications.

The structure of the paper is as follows: First, a short description of the investigated real test nanosized structure, its finite difference method approximation, and structure discretization are presented; then, the approximation scheme of the DPL model for the considered MEMS structure is proposed; after that, the Krylov subspace-based model order reduction technique is demonstrated; and finally, the simulation results obtained using both reduced and non-reduced versions of the DPL model are compared, analyzed, and discussed.

2. Mathematical Preliminaries of the Proposed Methodology

2.1. Finite Difference Method Approximation

Thermal analysis was carried out for the real microelectromechanical system (MEMS) nanostructure, presented in Figure 1, which was manufactured at the Polish Institute of Electron Technology [14,15]. This test structure includes two parallel platinum resistors, each 10 µm long. One resistor is treated as a heater, whereas the second resistor plays the role of a temperature sensor. The cross-sectional area of each resistor is 100 nm wide and 100 nm high. The distance between these resistors is 100 nm. The resistors are placed on a 100 nm wide silicon dioxide layer and both are stacked on a 500 µm thick silicon layer. A detailed description of the investigated structure, as well as the measurement process are found in [14,15].

In the investigated cross-sectional area of the MEMS test structure, the Dirichlet boundary conditions are at the bottom, while the Neumann boundary conditions are used on the left, right, and the top parts of the structure and their environment. The boundary conditions can be described by the following equations:

$$\begin{gathered} T_k\left(t\right)=0,t\in R_{+}{{\displaystyle \cup }}^{}\left\{0\right\}, \\ k\in \left\{1,2,3,\dots ,n_x\right\} \end{gathered}$$

(3)

$$\begin{gathered} {\mathit{q}}_k\left(t\right)=0,t\in R_{+}{{\displaystyle \cup }}^{}\left\{0\right\}, \\ k\in \left\{n_x+1,2\cdot n_x+1,\dots ,\left(n_y-1\right)\cdot n_x+1\right\}{{\displaystyle \cup }}^{}\left\{n_x,2\cdot n_x,\dots ,\left(n_y-1\right)\cdot n_x\right\}{{\displaystyle \cup }}^{} \\ {{\displaystyle \cup }}^{}\left\{\left(n_y-1\right)\cdot n_x+1,\left(n_y-1\right)\cdot n_x+2,\dots ,n_y\cdot n_x\right\} \end{gathered}$$

(4)

where n_x and n_y reflects the number of discretization nodes in both axes OX and OY, respectively.

Thermal simulations for the investigated MEMS structure were performed for its two-dimensional (2D) cross-sectional area, in the middle of the resistor. The thermal analysis of the structure was carried out using the DPL model. In order to make the analysis easier, the DPL model described by the system Equation (2) was equivalently transformed to the following 2D form [16]:

$$\begin{array}{c}{c}_v\cdot {\tau }_q\cdot \frac{{\partial }^{2}T\left(x,y,t\right)}{\partial t^2}+c_v\cdot \frac{\partial T\left(x,y,t\right)}{\partial t}-k\cdot {\tau }_T\cdot \frac{\partial \Delta T\left(x,y,t\right)}{\partial t}-k\cdot \Delta T\left(x,y,t\right)=\hfill \\ =q_V\left(x,y,t\right)+{\tau }_q\cdot \frac{\partial q_V\left(x,y,t\right)}{\partial t}\mathrm{for}x,y\in R,t\in R_{+}{{\displaystyle \cup }}^{}\left\{0\right\} \hfill \end{array}$$

(5)

As presented in Appendix A, the term ${\tau }_q\frac{\partial q_v\left(t\right)}{\partial t}$ can be neglected with an error of about −1.91%, for the simulation time t ≥ 3τ_T and $\left|err\right|\le 0.58\%$ for $t\ge 10{\tau }_T$. The $q_v\left(t+{\tau }_q\right)$ can also be used instead of $q_v\left(t\right)+{\tau }_q\frac{\partial q_v\left(t\right)}{\partial t}$ for $q_v\left(t\right)=c_1\cdot H\left(t\right)$, where H(t) is a Heaviside function, $c_1$ is a constant where $\left|c_1\right|<+\infty$. It is worthwhile highlighting that the presented values were obtained considering the temperature and heat flux time lags for platinum. The Laplacian ΔT was approximated using the finite difference method (FDM) approach for a rectangular mesh with a constant distance between nodes in both dimensions. The considered approximation is presented below:

$$\begin{array}{c}\Delta T\left(x,y,t\right)\approx \frac{T\left(x+\Delta x,y,t\right)+T\left(x,y+\Delta x,t\right)-4\cdot T\left(x,y,t\right)+T\left(x-\Delta x,y,t\right)+T\left(x,y-\Delta x,t\right)}{{\left(\Delta x\right)}^2},\hfill \\ \mathrm{dla}x,y\in R,t\in R_{+}{{\displaystyle \cup }}^{}\left\{0\right\} \hfill \end{array}$$

(6)

Then, the DPL Equation (5) becomes an ordinary differential equation (ODE) of a time variable. Such a constructed system of equations, including each investigated mesh node, was solved based on the backward differentiation formulas (BDF) approach [17,18,19,20], with initial conditions T(x,y,t) = 0 and ∂T(x,y,t)/∂t = 0 for t ≤ 0, which means zero initial conditions for Equation (16).

2.2. Description of Structure Discretization

The analyzed cross-sectional area of the MEMS structure was discretized using the mesh that can be described by the following equations [16]:

$${\mathit{q}}_k\left(t\right)=\mathit{q}\left(x,y,t\right),x=i\cdot \Delta x,y=j\cdot \Delta y$$

(7)

$$T_k\left(t\right)=T\left(x,y,t\right),x=i\cdot \Delta x,y=j\cdot \Delta y$$

(8)

i ∊ {1, 2, …, n_x}, j ∊ {1, 2, …, n_y}, k ∊ {1, 2, …, n_x ∙n_y}

In the Equations above, n_x and n_y reflect the number of discretization nodes in both axes. It means that the product n_x × n_y describes the number of nodes used in discretization mesh. Nodes are numbered from the left side to the right side. First, the first bottom row is considered. When all nodes in this row have already had their numbers, the procedure is repeated in the second row from the bottom side. This process is continued until all nodes in a top row are numbered. The graphical interpretation of the applied discretization mesh for the investigated cross-section of the MEMS structure, as well as nodes numbering approach, are demonstrated in Figure 2.

2.3. DPL Model Approximation Scheme for the Investigated MEMS Structure

Considering the imposed FDM assumptions, the DPL model can be expressed in the following way [22]:

$$\{\begin{array}{cc}\begin{array}{c}\mathit{M}\cdot \ddot{\mathit{T}}\left(t\right)+\mathit{D}\cdot \dot{\mathit{T}}\left(t\right)+\mathit{K}\cdot \mathit{T}\left(t\right)=\mathit{b}\cdot u\left(t\right)\\ \mathit{y}\left(t\right)={\mathit{c}}^T\cdot \mathit{T}\left(t\right)\end{array}& t\in R_{+}{{\displaystyle \cup }}^{}\left\{0\right\}\end{array}$$

(9)

where superscript T means a transposition. Moreover, T is the temperature function, while $\dot{\mathit{T}}$ and $\ddot{\mathit{T}}$ are its first and second time derivative, respectively. All of them are n_x ∙ n_y × 1 vectors. Other vectors and matrices can be described as follows [23]:

$$\begin{array}{c}\mathit{I},{\mathit{M}}_{\mathit{F}\mathit{D}\mathit{M}},\mathit{M},\mathit{D},\mathit{K},{\mathit{c}}^T\in R^{n_x\cdot n_y\times n_x\cdot n_y},\hfill \\ {\mathit{c}}_{\mathit{v}},\mathit{k},{\mathbf{\tau }}_{\mathit{q}},{\mathbf{\tau }}_{\mathit{T}},\mathit{b},\mathit{y}\in R^{n_x\cdot n_y\times 1},\hfill \\ u\in R\hfill \end{array}$$

$$\mathit{M}=diag\left({\mathbf{\tau }}_{\mathit{q}}\right) \circ diag\left({\mathit{c}}_{\mathit{v}}\right)\circ \mathit{I}$$

(10)

$$\mathit{D}=diag\left({\mathit{c}}_{\mathit{v}}\right)\circ \mathit{I}-\frac{1}{{\left(\Delta x\right)}^2}\cdot repmat\left({\mathbf{\tau }}_{\mathit{T}}\right)\circ repmat\left(\mathit{k}\right)\circ {\mathit{M}}_{\mathit{F}\mathit{D}\mathit{M}}$$

(11)

$$\mathit{K}=-\frac{1}{{\left(\Delta x\right)}^2}repmat\left(\mathit{k}\right)\circ {\mathit{M}}_{\mathit{F}\mathit{D}\mathit{M}}$$

(12)

$$\mathit{b}=a\cdot {\left[\begin{array}{ccccccccc}0& \cdots & 0& 1& \cdots & 1& 0& \cdots & 0\end{array}\right]}^T,a\in R_{+}$$

(13)

$$\mathit{c}=\mathit{I}$$

(14)

where I is identity matrix, an operator ◦ means a Hadamard product, function diag(∙) generates a diagonal matrix based on a given vector, and function repmat(∙) replicates a considered vector and generates a matrix of desired dimensions. In addition, non-zero elements of vector b indicate nodes with internal heat generation. Moreover, M_FDM is a matrix including FDM coefficients wgugc us determined based on Equation (6) and the considered boundary conditions. It is also worthwhile to note that the M_FDM coefficients describing the air area between platinum resistors’ surfaces, as well as the contact between them and the oxide, have been calculated considering a fractional order of the temperature function space derivative based on the Grünwald–Letnikov theory [22,24]. Thus, in this particular case, each one-dimensional part of Equation (6) is replaced by the following Equation [22]:

$$\begin{array}{c}\frac{1}{{\left(\Delta s\right)}^{{\alpha }_x}}\cdot {\displaystyle {\displaystyle \sum }_{k=0}^2}{\left(-1\right)}^k\frac{\Gamma \left({\alpha }_x+1\right)}{\Gamma \left(k+1\right)\cdot \Gamma \left({\alpha }_x-k+1\right)}T\left(\mathit{x}-k\cdot \Delta \mathit{x}+\frac{{\alpha }_x\cdot \Delta \mathit{x}}{2},t\right)\cdot 1=\hfill \\ =\frac{\left(\frac{{\alpha }_x}{2}-1\right)\cdot T\left(x+2\cdot \Delta x,y,t\right)+\left(2-\frac{{\alpha }_x}{2}-{\alpha }_x\cdot \left(\frac{{\alpha }_x}{2}-1\right)\right)\cdot T\left(x+\Delta x,y,t\right)}{{\left(\Delta x\right)}^{{\alpha }_x}}+\hfill \\ +\frac{\left(\frac{{\alpha }_x\cdot \left({\alpha }_x-1\right)}{2}\cdot \left(\frac{{\alpha }_x}{2}-1\right)-{\alpha }_x\cdot \left(2-\frac{{\alpha }_x}{2}\right)\right)\cdot T\left(x,y,t\right)+\left(\left(2-\frac{{\alpha }_x}{2}\right)\cdot \frac{{\alpha }_x\cdot \left({\alpha }_x-1\right)}{2}\right)\cdot T\left(x-\Delta x,y,t\right)}{{\left(\Delta x\right)}^{{\alpha }_x}}\hfill \\ +\frac{\left(\frac{{\alpha }_x}{2}-1\right)\cdot T\left(x,y+2\cdot \Delta x,t\right)+\left(2-\frac{{\alpha }_x}{2}-{\alpha }_x\cdot \left(\frac{{\alpha }_x}{2}-1\right)\right)\cdot T\left(x,y+\Delta x,t\right)}{{\left(\Delta x\right)}^{{\alpha }_x}}+\hfill \\ +\frac{\left(\frac{{\alpha }_x\cdot \left({\alpha }_x-1\right)}{2}\cdot \left(\frac{{\alpha }_x}{2}-1\right)-{\alpha }_x\cdot \left(2-\frac{{\alpha }_x}{2}\right)\right)\cdot T\left(x,y,t\right)+\left(\left(2-\frac{{\alpha }_x}{2}\right)\cdot \frac{{\alpha }_x\cdot \left({\alpha }_x-1\right)}{2}\right)\cdot T\left(x,y-\Delta x,t\right)}{{\left(\Delta x\right)}^{{\alpha }_x}}\hfill \\ \mathrm{for}{\alpha }_x\in \left(2,2.5\right),\Delta x\to 0 \hfill \end{array}$$

(15)

In order to make the analysis more clear, the system Equations (9) of the second-order equations has been transformed into the first-order Equations [21]:

$$\{\begin{array}{cc}\begin{array}{c}\mathit{E}\cdot \stackrel{·}{\overline{\mathit{T}}}\left(t\right)=\mathit{A}\cdot \overline{\mathit{T}}\left(t\right)+\mathit{B}\cdot u\left(t\right)\\ \mathit{y}\left(t\right)={\mathit{C}}^T\cdot \overline{\mathit{T}}\left(t\right)\end{array}& t\in R_{+}{{\displaystyle \cup }}^{}\left\{0\right\}\end{array}$$

(16)

where $\overline{\mathit{T}}$ and $\stackrel{·}{\overline{\mathit{T}}}$ are 2 ∙ n_x ∙ n_y × 1 vectors. First n_x ∙ n_y elements of $\overline{\mathit{T}}$ are the same, similar to the case of T vector, while its second part coincides with $\stackrel{·}{\mathit{T}}$. The first half of $\stackrel{·}{\overline{\mathit{T}}}$ includes elements of $\stackrel{·}{\mathit{T}}$, whereas its second half states $\ddot{\mathit{T}}$. In addition, the remaining matrices and vectors visible in Equation (16) present as follows [21,23,25]:

$$\mathit{E}=\left[\begin{array}{cc}\mathit{I}& \Theta \\ \Theta & \mathit{M}\end{array}\right],\mathit{A}=\left[\begin{array}{cc}\Theta & \mathit{I}\\ -\mathit{K}& -\mathit{D}\end{array}\right],\mathit{B}=\left[\begin{array}{c}{\Theta }_1\\ \mathit{b}\end{array}\right],\mathit{C}=\left[\begin{array}{cc}\mathit{c}& \Theta \\ \Theta & \mathit{c}\end{array}\right]$$

(17)

where Θ and Θ₁ are null matrices. The final solution is determined using the BDF method with a variable order between 1 and 5.

2.4. Krylov Subspace Method

As seen in Figure 2, the thermal characteristic of an electronic structure can refer to the construction of a system consisting of a high number of differential equations. Moreover, the larger the system the more accurate the results. However, preparation of a numerical solution of such a system can be very time-consuming and can require significant computational power. Thus, a solution for this problem is needed. One idea to deal with it is an order reduction of a constructed system of equations.

The order reduction process can be based on a moment matching method [26]. This method helps with moments calculation, being a negative coefficients of a system’s transfer function F_T in a Taylor’s series, around point 0 [26,27]. Its form can be as follows [25]:

$$F_T\left(l\right)=-{\displaystyle \sum }_{i=0}^{\infty }{m}_i\cdot l^i$$

(18)

where m_i means an ith moment. For example, for the system Equation (16), moments are determined according to the following Equation:

$$m_i={\mathit{C}}^T{\left({\mathit{A}}^{-1}\cdot \mathit{E}\right)}^i\cdot {\mathit{A}}^{-1}\cdot \mathit{B},i\in N{{\displaystyle \cup }}^{}\left\{0\right\}$$

(19)

The main idea of the model order reduction is related to finding a new system of equations, which is characterized by a significantly lower order than the original one. In this process, it is extremely important to consider the same transfer function F_T for both systems. It assures the existence of the same initial moments in the original and reduced systems of equations. However, due to the form of Equation (18), a direct numerical determination of such a proposed solution is impossible. Thus, to resolve this problem, the Krylov subspace method can be used.

In linear algebra, the Krylov space (or subspace) K of an order r for a certain n × n square matrix U and n-element vector s is a linear subspace of a space Rⁿ that is generated by the following vectors: s, U∙s, U²∙s, …, Uⁿ⁻¹∙s. Thus, the following Equation is true [11,28,29]:

$$K_r\left(\mathit{U},\mathit{s}\right)=\mathrm{span}\left\{\mathit{s},\mathit{U}\cdot \mathit{s},{\mathit{U}}^2\cdot \mathit{s},\dots ,{\mathit{U}}^{n-1}\cdot \mathit{s}\right\}$$

(20)

In the case of the previously analyzed system Equation (16), the matrix U is the same as A, while vector s is similar to B.

A Krylov subspace-based method is a numerical method that can find eigenvalues of large sparse matrices or solves a system with a high number of linear equations based on multiplications of vectors and matrices and operates on the determined vectors without the necessity of making additional operations on many matrices at the same time. Thus, starting with the vectors, the following vectors are determined, respectively: U∙s, U²∙s, etc.

Taking into consideration that consecutively determined vectors become linearly dependent relatively quickly, Krylov subspace-based methods require an additional application of some orthogonalization schemes. One of the most common of them are algorithms created by Lanczos [30,31,32,33,34], and Arnoldi [35,36,37,38,39,40,41]. The second algorithm mentioned was used in considerations presented in this paper. It can generate a set of orthonormal vectors which are determined using the modified Gram–Schmidt process (MGS) [42,43]. These vectors form a base Equation (20) of a Krylov subspace. Moreover, each vector from this base, starting with b, states consecutive columns of so-called transfer matrices V and W, which are used to determine coefficients matrices for reduced systems of linear equations. The algorithm stops after generating the first zero vectors, i.e., if the aggregated sum of absolute values of all vector’s coordinates does not exceed the tolerance value equal to 0.0001. The number of algorithm’s iterations (or the number of non-zero vectors) is equal to the reduced model order r. It is worthwhile highlighting that these zero vectors are not included in constructed V and W matrices.

In the investigated case, the V and W matrices are identical. Moreover, the number of iterations of the considered algorithm determines an order of the Krylov subspace and, at the same time, an order of a newly generated, reduced system of equations. For example, the reduced version of the system of Equation (16) is as follows [25]:

$$\{\begin{array}{cc}\begin{array}{c}{\mathit{E}}_{\mathit{r}}\cdot {\stackrel{·}{\overline{\mathit{T}}}}_{\mathit{r}}\left(t\right)={\mathit{A}}_{\mathit{r}}\cdot {\overline{\mathit{T}}}_{\mathit{r}}\left(t\right)+{\mathit{B}}_{\mathit{r}}\cdot u\left(t\right)\\ \mathit{y}\left(t\right)={\mathit{C}}_{\mathit{r}}\cdot {\overline{\mathit{T}}}_{\mathit{r}}\left(t\right)\end{array}& t\in R_{+}{{\displaystyle \cup }}^{}\left\{0\right\}\end{array}$$

(21)

where the system Equation (21) is solved using backward differentiation formulas (BDFs) of variable order between 1 and 5. Previous research has shown that it is one of the most effective methods for solving these types of equation systems 3.

3. Thermal Simulation Results

3.1. Pre-Simulation Assumptions

Thermal simulation of the investigated MEMS structure was carried out using Matlab environment and proposed approximation scheme for the DPL model. The simulation results were calculated using the following material thermal parameters [23,44] included in

Table 1. Considered material parameters’ values [23,44].

Material Name	$\mathit{k}\left[\frac{\mathit{W}}{\mathit{m}\cdot \mathit{K}}\right]$	$\mathit{\rho }\left[\frac{\mathbf{kg}}{{\mathit{m}}^3}\right]$	${\mathit{c}}_{\mathit{p}}\left[\frac{\mathit{J}}{\mathbf{kg}\cdot \mathit{K}}\right]$
Silicon (Si)	148	2 330	712
Silicon dioxide (SiO2)	1.38	2220	745
Platinum (Pt)	71.6	21 450	133
Platinum (Pt)	71.6	21 450	133
Air	0.0263	1.1614	1.007

:

Moreover, based on [23,45], the DPL model parameters were estimated. The platinum resistors are characterized by the heat flux time lag approximately equal to 550 ps and the temperature time lag established at 15.6 ns. For other materials, these values are set to 18 ns and 480 ns, respectively. In order to make an analysis easier, all results were normalized according to the following Equation:

$$T_k^{norm}\left(t\right)=\frac{T_k\left(t\right)}{\underset{t,k}{max}\left\{T_k\left(t\right)\right\}}k\in \left\{1,2,\dots ,n_x\cdot n_y\right\},t\in R_{+}{{\displaystyle \cup }}^{}\left\{0\right\}$$

(22)

3.2. Simulation Results

The thermal analyses of the investigated MEMS structure were divided into different areas. The first one, presented in this subsection, is related to the comparison between the results obtained using the reduced and non-reduced DPL model. In this area, the dynamic of average temperature rises in the platinum heater and the temperature sensor were investigated. Moreover, the temperature distribution in an entire cross-sectional area was considered for selected time points. It is also worthwhile highlighting that all results included in this subsection were received using a 10 nm distance between mesh nodes. Additionally, a comparison with FK and the measurement results was also carried out.

A comparison of the normalized average temperature rises over the time in the platinum heater and temperature sensor obtained using the reduced and non-reduced DPL model, the FK model, as well as real measurements is shown in Figure 3. It can be seen that there are almost no differences between the results plotted based on the reduced and non-reduced DPL models. The black solid line, which indicates the non-reduced DPL model results for the heater, coincide almost exactly with the red dashed line, which shows the outputs of the reduced DPL approach. A similar situation is also visible for the case of the temperature sensor. The black dotted line shows the results yielded using the non-reduced DPL model, and it coincides with the dashed blue curve, which shows outputs produced by the reduced DPL approach.

For comparison purposes, the measurement results which are indicated by the green lines were also plotted. It can be seen that the volatility over the time is very similar to the simulation results obtained using the DPL model, which confirms the correctness of the proposed approach.

The FK model produces significantly different outputs. The dashed and dotted black curve, which shows the temperature rise in the heater, as well as the dashed black line, which indicates the temperature changes over time in the temperature sensor, does not coincide with the DPL model or the measurement results. Thus, the FK model should not be used for temperature distribution determination at nanoscale.

In addition, an analysis of a temperature distribution inside an entire considered cross-sectional area of the MEMS structure was carried out. The results’ comparison, for selected time points, is demonstrated in Figure 4, Figure 5 and Figure 6.

For each investigated time point, i.e., at different stages of temperature rise, it can be seen that the differences between the results obtained using the reduced and non-reduced DPL model are almost unnoticeable, which also suggests a very good level of coincidence between both investigated DPL versions.

3.3. Relative Error Analysis

In order to assess the quality of results obtained using the reduced DPL model, an error analysis was carried out. Figure 7 shows a comparison of relative error values calculated between the reduced and non-reduced DPL model outputs for three selected points of time, similar to those investigated in the previous subsection. For all the mentioned time points, the relative error values in relation to each discretization mesh nodes were analyzed.

The maximum relative error value did not exceed a level of 5%, which confirms good accuracy of the reduced DPL model results. Moreover, the relative error decreased with an increase in the value of the investigated time points. For example, for higher observed temperature rise values, the relative error is at a level of 1%–1.35%. Furthermore, for steady state, a maximum relative error value is approximately equal to 0.01%. Such a small error value can be neglected. The decrease of the relative error value over time also suggests a convergence of the proposed approach.

3.4. Computational Complexity Analysis

Finally, a computational complexity of the non-reduced and reduced DPL models was investigated. First, an analysis of the time consumed during the solution calculation was carried out and is shown in Figure 8. It can be seen that the smaller the distance between mesh nodes the longer the time needed to obtain a temperature distribution in the investigated cross-sectional area of the MEMS structure. This is caused by a greater number of nodes in the smaller considered node distances.

The simulation time investigation in relation to the number of nodes in the analyzed cross-section is shown in Figure 9. The significant difference in simulation time needed for results preparation is visible. Considering the non-reduced DPL model, on the one hand, the simulation time can be estimated by a power function according to Equation (23). On the other hand, the time consumed using the reduced DPL model can be approximated based on a linear function presented in Equation (24). Moreover, statistics describing Equations for simulation time estimations are shown in

Table 2. Simulation time approximation.

Metrics	Adjusted R2	R2	RMSE
$t_{non-reduced}^{DPL}$	0.9988	0.9988	174.3
$t_{reduced}^{DPL}$	0.8299	0.829	0.06634

.

$$t_{non-reduced}^{DPL}\left(n\right)=1.96\cdot {10}^{-7}\cdot n^{2.544}$$

(23)

$$t_{reduced}^{DPL}\left(n\right)=2.546\cdot {10}^{-5}\cdot n-0.03782$$

(24)

where n reflects the number of nodes in the considered MEMS structure cross-section.

The values of metrics presented in Table 2 confirm the high quality of the prepared Equations for simulation time estimation. Moreover, based on Equations (23) and (24) and the simulation time analysis, it can be stated that the time complexity of an algorithm using the non-reduced DPL model is O(n^2.544), whereas the reduced DPL is characterized by a time complexity O(n). It is also worthwhile highlighting that the DPL system Equation (21) is linear (regarding the time variable), and therefore the Krylov subspace is calculated only once. Thus, it is not updated in each time step. Of course, the computational complexity of the method includes considerations regarding the Krylov subspace generation.

4. Conclusions

This paper includes an analysis of the quality and time complexity of two algorithms dedicated to solving heat transfer problems at nanoscale. The first one uses the modern DPL model which is a significant improvement as compared with one of the most common approaches based on the FK model. The second one also employs the DPL model, however in its reduced version. The DPL model order reduction is prepared based on the Krylov subspace method.

The analyses have shown that both the reduced and non-reduced DPL models produce high quality results that coincide with real measurements of the test structure. Moreover, the results obtained using the reduced DPL model are very similar to the results yielded based on the non-reduced DPL approach. It is also worthwhile highlighting that the relative error of approximation of temperature distribution determination inside the test structure, which was obtained using the reduced DPL model, was at a very low level. Furthermore, this error decreased over the time, which suggested a convergence of the proposed approach.

In addition, the reduced DPL model prepares a solution for a heat transfer problem significantly faster than the classic version of this heat transfer model. The time complexity of the non-reduced approach is O(n^2.544), whereas in the case of the reduced model, the complexity is O(n) only. Considering all the mentioned facts, it can be stated that the proposed approach obtained the high quality solution of the temperature distribution at nanoscale in a significantly shorter time than the classic approach, which is especially important to the future of designing and investigating advanced nanosized electronic structures.