INTEGRAL TRANSFORM SOLUTION OF HEAT
CONDUCTION IN ANISOTROPIC
HETEROGENEOUS MEDIA

Renato Cotta, Mohammed Lachi, Carolina Palma Naveira-Cotta, Arthur E.H.

Bruno

► To cite this version:

Renato Cotta, Mohammed Lachi, Carolina Palma Naveira-Cotta, Arthur E.H. Bruno. INTEGRAL TRANSFORM SOLUTION OF HEAT CONDUCTION IN ANISOTROPIC HETEROGENEOUS MEDIA. Computational Thermal Sciences: An International Journal, 2023, Computational Thermal Sciences: An International Journal, 16 (2), pp.65-90. 10.1615/ComputThermalScien.2023049260 . hal-04595064

HAL Id: hal-04595064

https://hal.science/hal-04595064

Submitted on 30 May 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
INTEGRAL TRANSFORM SOLUTION OF HEAT
CONDUCTION IN ANISOTROPIC
HETEROGENEOUS MEDIA

Renato Cotta, Mohammed Lachi, Carolina Palma Naveira-Cotta, Arthur E.H. Bruno

Bruno

► To cite this version:

Renato Cotta, Mohammed Lachi, Carolina Palma Naveira-Cotta, Arthur E.H. Bruno. INTEGRAL TRANSFORM SOLUTION OF HEAT CONDUCTION IN ANISOTROPIC HETEROGENEOUS MEDIA. Computational Thermal Sciences: An International Journal, 2023, Computational Thermal Sciences: An International Journal, 16 (2), pp.65-90. 10.1615/ComputThermalScien.2023049260 . hal-04595064

HAL Id: hal-04595064

https://hal.science/hal-04595064v1

Submitted on 30 May 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Integral Transform Solution of Heat Conduction in
Anisotropic Heterogeneous Media

Renato M. Cotta¹,², Mohammed Lachi³,
Carolina P. Naveira-Cotta¹, and Arthur E.H. Bruno¹

¹Laboratory of Nano and Microfluidics and Microsystems, LabMEMS, Mechanical Eng. Dept., POLI & COPPE, Federal University of Rio de Janeiro, UFRJ, RJ, Brazil

²Navy Research Institute, IPqM/CTMRJ, General Directorate of Nuclear and Technological Development, DGDNTM, Brazilian Navy, Rio de Janeiro, RJ, Brazil

³Matériaux et Ingénierie Mécanique (MATIM), Université de Reims Champagne-
Ardenne, CEDEX 2, 51687 Reims, France

*Corresponding author: cotta@mecanica.coppe.ufrj.br
Integral Transform Solution of Heat Conduction in
Anisotropic Heterogeneous Media

Renato M. Cotta¹,², Mohammed Lachi³,
Carolina P. Naveira-Cotta¹, and Arthur E.H. Bruno¹

¹Laboratory of Nano and Microfluidics and Microsystems, LabMEMS, Mechanical Eng. Dept., POLI & COPPE, Federal University of Rio de Janeiro, UFRJ, RJ, Brazil

²Navy Research Institute, IPqM/CTMRJ, General Directorate of Nuclear and Technological Development, DGDNTM, Brazilian Navy, Rio de Janeiro, RJ, Brazil

³Matériaux et Ingénierie Mécanique (MATIM), Université de Reims Champagne-
Ardenne, CEDEX 2, 51687 Reims, France

*Corresponding author: cotta@mecanica.coppe.ufrj.br
ABSTRACT

The Generalized Integral Transform Technique (GITT) is advanced to deal with conduction heat transfer in anisotropic heterogeneous media. A formal solution for exact integral transformation of conduction in anisotropic media is extended to account for heterogeneities expressed as space variable equation coefficients and source terms. The proposed eigenfunction expansion is based on biorthogonal eigenvalue problems which results in an exact integral transformation for linear problems and in a coupled transformed system for nonlinear situations. The associated eigenvalue problems are also handled through the GITT by considering simpler auxiliary eigenvalue problems of known analytical solution, leading to transformed algebraic eigenvalue problems. Finally, a two-dimensional transient test case is considered that presents an abrupt transition between isotropic and anisotropic materials yielding a marked change in thermal behavior in a defined region of interest formed by the anisotropic inclusion.

KEYWORDS

Conduction; heterogeneous media; anisotropic media; integral transforms; Generalized
Integral Transform Technique; eigenvalue problem; FGM; metamaterials
ABSTRACT

The Generalized Integral Transform Technique (GITT) is advanced to deal with conduction heat transfer in anisotropic heterogeneous media. A formal solution for exact integral transformation of conduction in anisotropic media is extended to account for heterogeneities expressed as space variable equation coefficients and source terms. The proposed eigenfunction expansion is based on biorthogonal eigenvalue problems which results in an exact integral transformation for linear problems and in a coupled transformed system for nonlinear situations. The associated eigenvalue problems are also handled through the GITT by considering simpler auxiliary eigenvalue problems of known analytical solution, leading to transformed algebraic eigenvalue problems. Finally, a two-dimensional transient test case is considered that presents an abrupt transition between isotropic and anisotropic materials yielding a marked change in thermal behavior in a defined region of interest formed by the anisotropic inclusion.

KEYWORDS

Conduction; heterogeneous media; anisotropic media; integral transforms; Generalized
Integral Transform Technique; eigenvalue problem; FGM; metamaterials
NOMENCLATURE

a,b Dimensions of test case rectangular substrate ($a = x_1 - x_0, b = y_1 - y_0$)

A, C, E Matrices in algebraic eigenvalue problem, Eq.(13a)

B Boundary conditions operator, Eq.(1e)

d Linear dissipation term coefficient, Eq.(1d)

f Initial condition function, Eq. (1b)

$\bar{f}_n$ Transformed initial condition, Eq. (6c)

$\bar{g}_n$ Transformed source term, Eq. (6d)

k Diffusion term coefficient, Eq. (1d)

L Space differential operator, Eq.(1d)

$M_k$ Normalization integral of auxiliary eigenvalue problem, Eq.(10b)

n Outward-drawn normal to the surface S, Eq. (1e)

$N_n$ Normalization integral (or norm), Eq. (4b)

P Nonlinear source term, Eq. (1a)

S Domain surface, Eq. (1c)

t Time variable or corresponding spatial variable, problem (1)

T Temperature or other potential (dependent variable), problem (1)

$\bar{T}_n$ Transformed temperature or potential, Eq.(5a)

$T_{in}$ Initial temperature in test-case

v Domain region, problem (1)

w Transient or convective term coefficient, Eq. (1a)

x Position vector, problem (1)
## NOMENCLATURE

a,b Dimensions of test case rectangular substrate ($a = x_1 - x_0, b = y_1 - y_0$)

A, C, E Matrices in algebraic eigenvalue problem, Eq.(13a)

B Boundary conditions operator, Eq.(1e)

d Linear dissipation term coefficient, Eq.(1d)

f Initial condition function, Eq. (1b)

$\bar{f}_n$ Transformed initial condition, Eq. (6c)

$\bar{g}_n$ Transformed source term, Eq. (6d)

k Diffusion term coefficient, Eq. (1d)

L Space differential operator, Eq.(1d)

$M_k$ Normalization integral of auxiliary eigenvalue problem, Eq.(10b)

n Outward-drawn normal to the surface S, Eq. (1e)

$N_n$ Normalization integral (or norm), Eq. (4b)

P Nonlinear source term, Eq. (1a)

S Domain surface, Eq. (1c)

t Time variable or corresponding spatial variable, problem (1)

T Temperature or other potential (dependent variable), problem (1)

$\bar{T}_n$ Transformed temperature or potential, Eq.(5a)

$T_{in}$ Initial temperature in test-case

V Domain region, problem (1)

w Transient or convective term coefficient, Eq. (1a)

x Position vector, problem (1)
$x,y$ Space coordinates in test-case, problem (19)

## Greek Symbols

$\alpha$ Boundary condition coefficient, Eq.(1e)

$\beta$ Boundary condition coefficient, Eq.(1e)

$\phi$ Nonlinear boundary source term defined in Eq. (1e)

$\mu_n$ Eigenvalues, problem (2)

$\psi_n$ Eigenfunctions, problem (2)

$\delta_{n,m}$ Kronecker delta operator, Eq. (4b)

$\lambda_k$ Auxiliary eigenvalues, problem (6)

$\Omega_k$ Auxiliary eigenfunctions, problem (6)

## Subscripts and Superscripts

$i,j$ Indices for tensorial notation, Eqs.(1d,e)

$n$ Index for eigenfunction/eigenvalue

$k,l$ Indices for auxiliary eigenfunction/eigenvalue

* Related to adjoint eigenvalue problem

^ Related to auxiliary eigenvalue problem

- Related to integral transformed function

~ Related to normalized eigenfunctions
$x,y$ Space coordinates in test-case, problem (19)

## Greek Symbols

$\alpha$ Boundary condition coefficient, Eq.(1e)

$\beta$ Boundary condition coefficient, Eq.(1e)

$\phi$ Nonlinear boundary source term defined in Eq. (1e)

$\mu_n$ Eigenvalues, problem (2)

$\psi_n$ Eigenfunctions, problem (2)

$\delta_{n,m}$ Kronecker delta operator, Eq. (4b)

$\lambda_k$ Auxiliary eigenvalues, problem (6)

$\Omega_k$ Auxiliary eigenfunctions, problem (6)

## Subscripts and Superscripts

$i,j$ Indices for tensorial notation, Eqs.(1d,e)

$n$ Index for eigenfunction/eigenvalue

$k,l$ Indices for auxiliary eigenfunction/eigenvalue

* Related to adjoint eigenvalue problem

^ Related to auxiliary eigenvalue problem

- Related to integral transformed function

~ Related to normalized eigenfunctions
# INTRODUCTION

Heat conduction in anisotropic media is directly associated with applications dealing with composite materials, laminated metals, wood, crystals and several other natural or man-made materials [1-4]. Often, such materials are inherently heterogeneous or are assembled together with isotropic substrates or connecting elements that participate in the heat transfer process, which in either case can be represented as spatially variable thermophysical properties, accounting for anisotropy wherever required. Such is the case, for instance, in a class of anisotropic functionally graded materials (FGM) [5-6] which present a tailored marked spatial variation of the thermophysical properties, aimed at achieving a specific functionality. Also recently, a class of thermal metamaterials has been advanced [7-10], which are designed based on an extended optical transformation theory that results in an anisotropic embedded region for achieving a desired controlled thermal behavior in a certain sub-domain, such as cloaking, rotation, concentration, camouflage, or illusion.

Exact solutions for transient heat conduction with the tensorial representation of the thermal conductivity are readily available for linear special cases such as in infinite or semi-infinite media or homogeneous orthotropic materials [1], but less evident for general anisotropy, heterogeneity, and finite multidimensional media. A generalization of the classical integral transform method [11] was proposed in [12], which accounts for general non-symmetric anisotropy in finite regions. Since through the classical integral transformation procedure [11], the resulting eigenvalue problem would be non-self adjoint for this more general situation, an adjoint eigenvalue problem is proposed in [12] that forms a biorthogonal set with the original one. Therefore, both the anisotropy and heterogeneity effects are transferred to the adjoint eigenvalue problems, which become fairly involved for the general case, though a formal exact solution is provided for the temperature distribution. This study was later on extended through a convenient matrix operator form [13], covering a wide class of linear diffusion problems.

Along the years and most markedly after the unification of the integral transform method in seven different classes of linear problems in heat and mass diffusion [11], a hybrid numerical-analytical extension to the classical approach, known as the Generalized Integral Transform Technique (GITT) [14] has been progressively advanced, as reviewed in different sources [15-24], aimed at handling different classes of formulations that would not allow for an exact integral transformation in the sense of [11].
# INTRODUCTION

Heat conduction in anisotropic media is directly associated with applications dealing with composite materials, laminated metals, wood, crystals and several other natural or man-made materials [1-4]. Often, such materials are inherently heterogeneous or are assembled together with isotropic substrates or connecting elements that participate in the heat transfer process, which in either case can be represented as spatially variable thermophysical properties, accounting for anisotropy wherever required. Such is the case, for instance, in a class of anisotropic functionally graded materials (FGM) [5-6] which present a tailored marked spatial variation of the thermophysical properties, aimed at achieving a specific functionality. Also recently, a class of thermal metamaterials has been advanced [7-10], which are designed based on an extended optical transformation theory that results in an anisotropic embedded region for achieving a desired controlled thermal behavior in a certain sub-domain, such as cloaking, rotation, concentration, camouflage, or illusion.

Exact solutions for transient heat conduction with the tensorial representation of the thermal conductivity are readily available for linear special cases such as in infinite or semi-infinite media or homogeneous orthotropic materials [1], but less evident for general anisotropy, heterogeneity, and finite multidimensional media. A generalization of the classical integral transform method [11] was proposed in [12], which accounts for general non-symmetric anisotropy in finite regions. Since through the classical integral transformation procedure [11], the resulting eigenvalue problem would be non-self adjoint for this more general situation, an adjoint eigenvalue problem is proposed in [12] that forms a biorthogonal set with the original one. Therefore, both the anisotropy and heterogeneity effects are transferred to the adjoint eigenvalue problems, which become fairly involved for the general case, though a formal exact solution is provided for the temperature distribution. This study was later on extended through a convenient matrix operator form [13], covering a wide class of linear diffusion problems.

Along the years and most markedly after the unification of the integral transform method in seven different classes of linear problems in heat and mass diffusion [11], a hybrid numerical-analytical extension to the classical approach, known as the Generalized Integral Transform Technique (GITT) [14] has been progressively advanced, as reviewed in different sources [15-24], aimed at handling different classes of formulations that would not allow for an exact integral transformation in the sense of [11].
Thus, diffusion and convection-diffusion problems with time-dependent coefficients, nonlinear problems, irregular geometries, non-classical eigenvalue problems, moving boundary problems, boundary layer equations, and Navier-Stokes equations, among others, were systematically approached and properly handled through the more flexible hybrid methodology. A few advancements allowed for by the generalized approach are particularly relevant to the present work, the first one being the solution of eigenvalue problems by the GITT itself, as proposed in [25], and its application in the analysis of diffusion in heterogeneous media [26-27]. This advancement was particularly important in allowing for the solution of fairly general eigenvalue problems with arbitrarily variable coefficients, by considering a simpler auxiliary eigenvalue problem as the expansion base, thus transforming cumbersome differential eigenvalue problems into readily solvable algebraic ones. A second advancement is related to a single domain reformulation strategy [28], which allows for the representation of multiregion problems as one single formulation with spatially variable coefficients accounting for the different materials transitions. This approach was then applied to conjugated conduction-convection problems [29-31] and heat transfer in heterogeneous and complex media [32-33], among other applications in heat and fluid flow.

The present work is concerned with revisiting the formal exact solution of [12] in dealing with heat conduction in heterogeneous anisotropic media, while extending it to nonlinear problems in the framework of the GITT, as first generalized in [19]. The solution of the biorthogonal eigenvalue problems is then proposed through the GITT itself, considering a pair of simpler biorthogonal auxiliary eigenvalue problem in the eigenfunctions expansion in solving the original more involved problems. For the sake of completeness, the traditional GITT formalism is also recalled, as an alternative solution path for heterogeneous anisotropic media, accounting for cross derivatives as source terms, and directly handling the integral transformation with an adjoint eigenvalue problem with space variable coefficients. A transient two-dimensional test case represented by an anisotropic material inclusion in an isotropic rectangular substrate is considered for illustrating the convergence behavior of the proposed integral transform solution. A well-established commercial finite element method software [34] is also considered for comparison purposes and the present error-controlled solutions provide benchmark results for its convergence analysis.
Thus, diffusion and convection-diffusion problems with time-dependent coefficients, nonlinear problems, irregular geometries, non-classical eigenvalue problems, moving boundary problems, boundary layer equations, and Navier-Stokes equations, among others, were systematically approached and properly handled through the more flexible hybrid methodology. A few advancements allowed for by the generalized approach are particularly relevant to the present work, the first one being the solution of eigenvalue problems by the GITT itself, as proposed in [25], and its application in the analysis of diffusion in heterogeneous media [26-27]. This advancement was particularly important in allowing for the solution of fairly general eigenvalue problems with arbitrarily variable coefficients, by considering a simpler auxiliary eigenvalue problem as the expansion base, thus transforming cumbersome differential eigenvalue problems into readily solvable algebraic ones. A second advancement is related to a single domain reformulation strategy [28], which allows for the representation of multiregion problems as one single formulation with spatially variable coefficients accounting for the different materials transitions. This approach was then applied to conjugated conduction-convection problems [29-31] and heat transfer in heterogeneous and complex media [32-33], among other applications in heat and fluid flow.

The present work is concerned with revisiting the formal exact solution of [12] in dealing with heat conduction in heterogeneous anisotropic media, while extending it to nonlinear problems in the framework of the GITT, as first generalized in [19]. The solution of the biorthogonal eigenvalue problems is then proposed through the GITT itself, considering a pair of simpler biorthogonal auxiliary eigenvalue problem in the eigenfunctions expansion in solving the original more involved problems. For the sake of completeness, the traditional GITT formalism is also recalled, as an alternative solution path for heterogeneous anisotropic media, accounting for cross derivatives as source terms, and directly handling the integral transformation with an adjoint eigenvalue problem with space variable coefficients. A transient two-dimensional test case represented by an anisotropic material inclusion in an isotropic rectangular substrate is considered for illustrating the convergence behavior of the proposed integral transform solution. A well-established commercial finite element method software [34] is also considered for comparison purposes and the present error-controlled solutions provide benchmark results for its convergence analysis.
# ANALYSIS

## Problem formulation

Transient heat conduction in a general anisotropic finite medium is considered, where heterogeneities are accounted for by arbitrarily space variable thermophysical properties and source terms, defined in region V with bounding surface S. The temperature distribution is governed by [12]:

$$w(\mathbf{x}) \frac{\partial T(\mathbf{x}, t)}{\partial t} + LT(\mathbf{x}, t) = P(\mathbf{x}, t, T), \quad \text{in } \mathbf{x} \in V, t > 0 \quad (1a)$$

with initial and boundary conditions given as

$$T(\mathbf{x}, 0) = f(\mathbf{x}), \quad \mathbf{x} \in V \quad (1b)$$

$$BT(\mathbf{x}, t) = \phi(\mathbf{x}, t, T), \quad \mathbf{x} \in S, t > 0 \quad (1c)$$

where the equation and boundary operators are written in tensor notation as

$$L = - \frac{\partial}{\partial x_i} \left( k_{ij}(\mathbf{x}) \frac{\partial}{\partial x_j} \right) + d(\mathbf{x}) \quad (1d,e)$$

$$B = \alpha(\mathbf{x}) + \beta(\mathbf{x}) n_i k_{ij}(\mathbf{x}) \frac{\partial}{\partial x_j}$$

where $\alpha(\mathbf{x})$ and $\beta(\mathbf{x})$ are prescribed coefficients that define the boundary conditions kinds and $n_i$ are the components of the outward drawn normal to boundary surface S in each coordinate direction $x_i$ ($i=1, 2, 3$).

It should be recalled that system (1) is in fact more general than it may seem, since any nonlinearities in the equation or boundary condition coefficients can be collected into the nonlinear source terms, $P(\mathbf{x}, t, T)$ and $\phi(\mathbf{x}, t, T)$, while characteristic $\mathbf{x}$-dependent coefficients are kept in the original operators. For isotropic media the system is simplified with $k_{ij}(\mathbf{x}) = k(\mathbf{x})$ and the formal solution to problem (1) is readily given in [14,16] through the Generalized Integral Transform Technique (GITT). For linear source terms, $P(\mathbf{x}, t)$ and $\phi(\mathbf{x}, t)$, and isotropic media, system (1) reduces to the so called Class I problems in [11] and formal exact solutions are readily available by the integral transform method.

## Biorthogonal eigenfunction expansion
# ANALYSIS

## Problem formulation

Transient heat conduction in a general anisotropic finite medium is considered, where heterogeneities are accounted for by arbitrarily space variable thermophysical properties and source terms, defined in region V with bounding surface S. The temperature distribution is governed by [12]:

$$w(\mathbf{x}) \frac{\partial T(\mathbf{x}, t)}{\partial t} + LT(\mathbf{x}, t) = P(\mathbf{x}, t, T), \quad \text{in } \mathbf{x} \in V, t > 0 \quad (1a)$$

with initial and boundary conditions given as

$$T(\mathbf{x}, 0) = f(\mathbf{x}), \quad \mathbf{x} \in V \quad (1b)$$

$$BT(\mathbf{x}, t) = \phi(\mathbf{x}, t, T), \quad \mathbf{x} \in S, t > 0 \quad (1c)$$

where the equation and boundary operators are written in tensor notation as

$$L = - \frac{\partial}{\partial x_i} \left( k_{ij}(\mathbf{x}) \frac{\partial}{\partial x_j} \right) + d(\mathbf{x}) \quad (1d,e)$$

$$B = \alpha(\mathbf{x}) + \beta(\mathbf{x}) n_i k_{ij}(\mathbf{x}) \frac{\partial}{\partial x_j}$$

where $\alpha(\mathbf{x})$ and $\beta(\mathbf{x})$ are prescribed coefficients that define the boundary conditions kinds and $n_i$ are the components of the outward drawn normal to boundary surface S in each coordinate direction $x_i$ ($i=1, 2, 3$).

It should be recalled that system (1) is in fact more general than it may seem, since any nonlinearities in the equation or boundary condition coefficients can be collected into the nonlinear source terms, $P(\mathbf{x}, t, T)$ and $\phi(\mathbf{x}, t, T)$, while characteristic $\mathbf{x}$-dependent coefficients are kept in the original operators. For isotropic media the system is simplified with $k_{ij}(\mathbf{x}) = k(\mathbf{x})$ and the formal solution to problem (1) is readily given in [14,16] through the Generalized Integral Transform Technique (GITT). For linear source terms, $P(\mathbf{x}, t)$ and $\phi(\mathbf{x}, t)$, and isotropic media, system (1) reduces to the so called Class I problems in [11] and formal exact solutions are readily available by the integral transform method.

## Biorthogonal eigenfunction expansion
Also for linear source terms, but for any general anisotropic thermal conductivity, a formal exact solution through the integral transform method has been proposed [12,13], based on a biorthogonal set of eigenfunctions. Here, we first briefly show how this formal solution can be extended to the nonlinear situation in eqs.(1). After dropping out the source terms in both the heat conduction equation and the boundary condition, the eigenvalue problem obtained from separation of variables is given by:

$$L\psi(\mathbf{x}) = \mu^2 w(\mathbf{x})\psi(\mathbf{x}), \quad \text{in } \mathbf{x} \in V \qquad (2a)$$

with boundary conditions given as

$$B\psi(\mathbf{x}) = 0, \quad \mathbf{x} \in S \qquad (2b)$$

However, for the general anisotropy tensor here considered, the eigenvalue problem given by eqs. (2) is non-self-adjoint and the resulting eigenfunctions are not orthogonal, thus the classical approach in [11,14] is not directly applicable. Then, following the formalism in [12-13,19], an adjoint eigenvalue problem to problem (2) is proposed, which will provide eigenfunctions that are orthogonal to those from problem (2), given by:

$$L^*\psi^*(\mathbf{x}) = \mu^2 w(\mathbf{x})\psi^*(\mathbf{x}), \quad \text{in } \mathbf{x} \in V \qquad (3a)$$

with boundary conditions given as

$$B^*\psi^*(\mathbf{x}) = 0, \quad \mathbf{x} \in S \qquad (3b)$$

where the adjoint operators to L and B are given by

$$L^* = -\frac{\partial}{\partial x_i} \left( k_{ji}(\mathbf{x}) \frac{\partial}{\partial x_j} \right) + d(\mathbf{x}) \qquad (3c,d)$$
$$B^* = \alpha(\mathbf{x}) + \beta(\mathbf{x}) n_i k_{ji}(\mathbf{x}) \frac{\partial}{\partial x_j}$$

It is shown in [12] that the eigenfunctions $\psi(\mathbf{x})$ and $\psi^*(\mathbf{x})$ obey the following orthogonality condition, forming a biorthogonal set:

$$\int_V w(\mathbf{x}) \psi_n(\mathbf{x}) \psi_m^*(\mathbf{x}) d\mathbf{v} = \delta_{n,m} N_n \qquad (4a)$$

where $\delta_{n,m}$ is the Kronecker delta and the norms are given by

$$N_n = \int_V w(\mathbf{x}) \psi_n(\mathbf{x}) \psi_n^*(\mathbf{x}) d\mathbf{v} \qquad (4b)$$

From this property, one may readily construct the integral transform-inverse pair, given by:
Also for linear source terms, but for any general anisotropic thermal conductivity, a formal exact solution through the integral transform method has been proposed [12,13], based on a biorthogonal set of eigenfunctions. Here, we first briefly show how this formal solution can be extended to the nonlinear situation in eqs.(1). After dropping out the source terms in both the heat conduction equation and the boundary condition, the eigenvalue problem obtained from separation of variables is given by:

$$L\psi(\mathbf{x}) = \mu^2 w(\mathbf{x})\psi(\mathbf{x}), \quad \text{in } \mathbf{x} \in V \qquad (2a)$$

with boundary conditions given as

$$B\psi(\mathbf{x}) = 0, \quad \mathbf{x} \in S \qquad (2b)$$

However, for the general anisotropy tensor here considered, the eigenvalue problem given by eqs. (2) is non-self-adjoint and the resulting eigenfunctions are not orthogonal, thus the classical approach in [11,14] is not directly applicable. Then, following the formalism in [12-13,19], an adjoint eigenvalue problem to problem (2) is proposed, which will provide eigenfunctions that are orthogonal to those from problem (2), given by:

$$L^*\psi^*(\mathbf{x}) = \mu^2 w(\mathbf{x})\psi^*(\mathbf{x}), \quad \text{in } \mathbf{x} \in V \qquad (3a)$$

with boundary conditions given as

$$B^*\psi^*(\mathbf{x}) = 0, \quad \mathbf{x} \in S \qquad (3b)$$

where the adjoint operators to L and B are given by

$$L^* = -\frac{\partial}{\partial x_i} \left( k_{ji}(\mathbf{x}) \frac{\partial}{\partial x_j} \right) + d(\mathbf{x}) \qquad (3c,d) \\ B^* = \alpha(\mathbf{x}) + \beta(\mathbf{x}) n_i k_{ji}(\mathbf{x}) \frac{\partial}{\partial x_j}$$

It is shown in [12] that the eigenfunctions $\psi(\mathbf{x})$ and $\psi^*(\mathbf{x})$ obey the following orthogonality condition, forming a biorthogonal set:

$$\int_V w(\mathbf{x})\psi_n(\mathbf{x})\psi_m^*(\mathbf{x})d\mathbf{v} = \delta_{n,m}N_n \qquad (4a)$$

where $\delta_{n,m}$ is the Kronecker delta and the norms are given by

$$N_n = \int_V w(\mathbf{x})\psi_n(\mathbf{x})\psi_n^*(\mathbf{x})d\mathbf{v} \qquad (4b)$$

From this property, one may readily construct the integral transform-inverse pair, given by:
$$ \bar{T}_n(t) = \int_V w(\mathbf{x}) T(\mathbf{x}, t) \psi_n^*(\mathbf{x}) dv \quad \text{transform} \qquad (5a) $$

$$ T(\mathbf{x}, t) = \sum_{n=1}^{\infty} \frac{1}{N_n} \psi_n(\mathbf{x}) \bar{T}_n(t) \quad \text{inverse} \qquad (5b) $$

Then, the integral transformation of eq.(1a) is performed employing the operator $\int_V (.) \psi_n^*(\mathbf{x}) dv$, following the formalism in [12-13], to yield the transformed system below:

$$ \frac{d\bar{T}_n(t)}{dt} + \mu_n^2 \bar{T}_n(t) = \bar{g}_n(t, \bar{T}(t)), \quad n = 1, 2, \dots, t > 0 \qquad (6a) $$

$$ \bar{T}_n(0) = \bar{f}_n \qquad (6b) $$

where the transformed source terms and initial conditions are given by

$$ \bar{f}_n = \int_V w(\mathbf{x}) f(\mathbf{x}) \psi_n^*(\mathbf{x}) dv \qquad (6c) $$

$$ \bar{g}_n(t, \bar{T}(t)) = \int_S \phi(\mathbf{x}, t, T) \left[ \frac{\psi_n^*(\mathbf{x}) - n_t k_{ji}(\mathbf{x}) \frac{\partial \psi_n^*(\mathbf{x})}{\partial x_j}}{\alpha(\mathbf{x}) + \beta(\mathbf{x})} \right] ds + \int_V P(\mathbf{x}, t, T) \psi_n^*(\mathbf{x}) dv \qquad (6d) $$

The transformed ordinary differential system (6), in the more general nonlinear case, can be solved through well tested routines for initial value problems, such as the function NDSolve in the Mathematica system [35]. For linear source terms, the solution of the transformed ordinary differential system (6) can be obtained analytically, as shown in [11,16], as:

$$ \bar{T}_n(t) = \bar{f}_n e^{-\mu_n^2 t} + \int_0^t \bar{g}_n(t') e^{-\mu_n^2 (t-t')} dt' \qquad (7) $$

In either case, the inverse formulae (5b) is then recalled to reconstruct the desired temperature distribution, a fully analytical solution for the linear case [12-13], or hybrid numerical-analytical for the nonlinear case.

## GITT solution of biorthogonal eigenvalue problems

The above formal solution requires the computation of eigenvalues and eigenfunctions for both problems (2) and (3), which might not be feasible in exact form. Therefore, the GITT approach itself can be employed in transforming the original differential eigenvalue problems into algebraic ones, by proposing a simpler auxiliary eigenvalue problem to provide a base for the eigenfunction expansions, such as described
$$ \bar{T}_n(t) = \int_V w(\mathbf{x})T(\mathbf{x}, t)\psi_n^*(\mathbf{x})dv \quad \text{transform} \qquad (5a) $$

$$ T(\mathbf{x}, t) = \sum_{n=1}^{\infty} \frac{1}{N_n} \psi_n(\mathbf{x}) \bar{T}_n(t) \quad \text{inverse} \qquad (5b) $$

Then, the integral transformation of eq.(1a) is performed employing the operator $\int_V (.) \psi_n^*(\mathbf{x})dv$, following the formalism in [12-13], to yield the transformed system below:

$$ \frac{d\bar{T}_n(t)}{dt} + \mu_n^2 \bar{T}_n(t) = \bar{g}_n(t, \bar{T}(t)), \quad n = 1, 2, \dots, t > 0 \qquad (6a) $$

$$ \bar{T}_n(0) = \bar{f}_n \qquad (6b) $$

where the transformed source terms and initial conditions are given by

$$ \bar{f}_n = \int_V w(\mathbf{x}) f(\mathbf{x}) \psi_n^*(\mathbf{x}) dv \qquad (6c) $$

$$ \bar{g}_n(t, \bar{T}(t)) = \int_S \phi(\mathbf{x}, t, T) \left[ \frac{\psi_n^*(\mathbf{x}) - n_t k_{ji}(\mathbf{x}) \frac{\partial \psi_n^*(\mathbf{x})}{\partial x_j}}{\alpha(\mathbf{x}) + \beta(\mathbf{x})} \right] ds + \int_V P(\mathbf{x}, t, T) \psi_n^*(\mathbf{x}) dv \qquad (6d) $$

The transformed ordinary differential system (6), in the more general nonlinear case, can be solved through well tested routines for initial value problems, such as the function NDSolve in the Mathematica system [35]. For linear source terms, the solution of the transformed ordinary differential system (6) can be obtained analytically, as shown in [11,16], as:

$$ \bar{T}_n(t) = \bar{f}_n e^{-\mu_n^2 t} + \int_0^t \bar{g}_n(t') e^{-\mu_n^2 (t-t')} dt' \qquad (7) $$

In either case, the inverse formulae (5b) is then recalled to reconstruct the desired temperature distribution, a fully analytical solution for the linear case [12-13], or hybrid numerical-analytical for the nonlinear case.

## GITT solution of biorthogonal eigenvalue problems

The above formal solution requires the computation of eigenvalues and eigenfunctions for both problems (2) and (3), which might not be feasible in exact form. Therefore, the GITT approach itself can be employed in transforming the original differential eigenvalue problems into algebraic ones, by proposing a simpler auxiliary eigenvalue problem to provide a base for the eigenfunction expansions, such as described
in [16,23,25]. According to the more traditional approach in solving eigenvalue problems, a simpler self-adjoint auxiliary problem could be chosen to provide the base of the eigenfunction expansion in solving problems (2) and (3). However, a similar path of considering a biorthogonal expansion is here followed also in solving the eigenvalue problems. In fact, the proposed solution that follows will automatically reduce to the traditional one once a self-adjoint eigenvalue problem is selected to be the auxiliary problem.

Thus, consider the following adjoint auxiliary eigenvalue problems with simpler coefficients, $\hat{w}(\mathbf{x})$, $\hat{k}_{ij}(\mathbf{x})$, $\hat{d}(\mathbf{x})$, $\hat{\alpha}(\mathbf{x})$, $\hat{\beta}(\mathbf{x})$, so as to allow for analytical solutions:

$$ \hat{L}\Omega(\mathbf{x}) = \lambda^2 \hat{w}(\mathbf{x})\Omega(\mathbf{x}), \quad \text{in } \mathbf{x} \in V \qquad (8a) $$

$$ \hat{B}\Omega(\mathbf{x}) = 0, \quad \mathbf{x} \in S \qquad (8b) $$

where the equation and boundary operators are written as

$$ \begin{aligned} \hat{L} &\equiv -\frac{\partial}{\partial x_i} \left( \hat{k}_{ij}(\mathbf{x}) \frac{\partial}{\partial x_j} \right) + \hat{d}(\mathbf{x}) \\ \hat{B} &\equiv \hat{\alpha}(\mathbf{x}) + \hat{\beta}(\mathbf{x}) n_i \hat{k}_{ij}(\mathbf{x}) \frac{\partial}{\partial x_j} \end{aligned} \qquad (8c,d) $$

and,

$$ \hat{L}^* \Omega^*(\mathbf{x}) = \lambda^2 \hat{w}(\mathbf{x}) \Omega^*(\mathbf{x}), \quad \text{in } \mathbf{x} \in V \qquad (9a) $$

$$ \hat{B}^* \Omega^*(\mathbf{x}) = 0, \quad \mathbf{x} \in S \qquad (9b) $$

where the adjoint operators to $\hat{L}$ and $\hat{B}$ are given by

$$ \begin{aligned} \hat{L}^* &\equiv -\frac{\partial}{\partial x_i} \left( \hat{k}_{ji}(\mathbf{x}) \frac{\partial}{\partial x_j} \right) + \hat{d}(\mathbf{x}) \\ \hat{B}^* &\equiv \hat{\alpha}(\mathbf{x}) + \hat{\beta}(\mathbf{x}) n_i \hat{k}_{ji}(\mathbf{x}) \frac{\partial}{\partial x_j} \end{aligned} \qquad (9c,d) $$

Thus, the auxiliary eigenfunctions $\Omega(\mathbf{x})$ and $\Omega^*(\mathbf{x})$ obey the following orthogonality condition, forming a biorthogonal set:

$$ \int_V \hat{w}(\mathbf{x}) \Omega_k(\mathbf{x}) \Omega_l^*(\mathbf{x}) d\mathbf{v} = \delta_{k,l} M_k \qquad (10a) $$

where $\delta_{k,l}$ is the Kronecker delta and the norms are given by

$$ M_k = \int_V \hat{w}(\mathbf{x}) \Omega_k(\mathbf{x}) \Omega_k^*(\mathbf{x}) d\mathbf{v} \qquad (10b) $$
in [16,23,25]. According to the more traditional approach in solving eigenvalue problems, a simpler self-adjoint auxiliary problem could be chosen to provide the base of the eigenfunction expansion in solving problems (2) and (3). However, a similar path of considering a biorthogonal expansion is here followed also in solving the eigenvalue problems. In fact, the proposed solution that follows will automatically reduce to the traditional one once a self-adjoint eigenvalue problem is selected to be the auxiliary problem.

Thus, consider the following adjoint auxiliary eigenvalue problems with simpler coefficients, $\hat{w}(\mathbf{x})$, $\hat{k}_{ij}(\mathbf{x})$, $\hat{d}(\mathbf{x})$, $\hat{\alpha}(\mathbf{x})$, $\hat{\beta}(\mathbf{x})$, so as to allow for analytical solutions:

$$ \hat{L}\Omega(\mathbf{x}) = \lambda^2 \hat{w}(\mathbf{x})\Omega(\mathbf{x}), \quad \text{in } \mathbf{x} \in V \qquad (8a) $$

$$ \hat{B}\Omega(\mathbf{x}) = 0, \quad \mathbf{x} \in S \qquad (8b) $$

where the equation and boundary operators are written as

$$ \begin{aligned} \hat{L} &\equiv -\frac{\partial}{\partial x_i} \left( \hat{k}_{ij}(\mathbf{x}) \frac{\partial}{\partial x_j} \right) + \hat{d}(\mathbf{x}) & (8c,d) \\ \hat{B} &\equiv \hat{\alpha}(\mathbf{x}) + \hat{\beta}(\mathbf{x}) n_i \hat{k}_{ij}(\mathbf{x}) \frac{\partial}{\partial x_j} \end{aligned} $$

and,

$$ \hat{L}^* \Omega^*(\mathbf{x}) = \lambda^2 \hat{w}(\mathbf{x}) \Omega^*(\mathbf{x}), \quad \text{in } \mathbf{x} \in V \qquad (9a) $$

$$ \hat{B}^* \Omega^*(\mathbf{x}) = 0, \quad \mathbf{x} \in S \qquad (9b) $$

where the adjoint operators to $\hat{L}$ and $\hat{B}$ are given by

$$ \begin{aligned} \hat{L}^* &\equiv -\frac{\partial}{\partial x_i} \left( \hat{k}_{ji}(\mathbf{x}) \frac{\partial}{\partial x_j} \right) + \hat{d}(\mathbf{x}) & (9c,d) \\ \hat{B}^* &\equiv \hat{\alpha}(\mathbf{x}) + \hat{\beta}(\mathbf{x}) n_i \hat{k}_{ji}(\mathbf{x}) \frac{\partial}{\partial x_j} \end{aligned} $$

Thus, the auxiliary eigenfunctions $\Omega(\mathbf{x})$ and $\Omega^*(\mathbf{x})$ obey the following orthogonality condition, forming a biorthogonal set:

$$ \int_V \hat{w}(\mathbf{x}) \Omega_k(\mathbf{x}) \Omega_l^*(\mathbf{x}) d\mathbf{v} = \delta_{k,l} M_k \qquad (10a) $$

where $\delta_{k,l}$ is the Kronecker delta and the norms are given by

$$ M_k = \int_V \hat{w}(\mathbf{x}) \Omega_k(\mathbf{x}) \Omega_k^*(\mathbf{x}) d\mathbf{v} \qquad (10b) $$
Then, the transform-inverse pair is proposed for the desired eigenfunction, $\psi_n(\mathbf{x})$

as:

$$ \bar{\psi}_{n,k} = \int_V \hat{w}(\mathbf{x}) \Omega_k^*(\mathbf{x}) \psi_n(\mathbf{x}) dv \quad \text{transform} \qquad (11a) $$

$$ \psi_n(\mathbf{x}) = \sum_{k=1}^{\infty} \frac{1}{M_k} \Omega_k(\mathbf{x}) \bar{\psi}_{n,k} \quad \text{inverse} \qquad (11b) $$

Similarly, the transform-inverse pair for the adjoint eigenfunction, $\psi_n^*(\mathbf{x})$ becomes:

$$ \bar{\psi}_{n,k}^* = \int_V \hat{w}(\mathbf{x}) \Omega_k(\mathbf{x}) \psi_n^*(\mathbf{x}) dv \quad \text{transform} \qquad (12a) $$

$$ \psi_n^*(\mathbf{x}) = \sum_{k=1}^{\infty} \frac{1}{M_k} \Omega_k^*(\mathbf{x}) \bar{\psi}_{n,k}^* \quad \text{inverse} \qquad (12b) $$

Now, we proceed with the integral transformation of eq.(2a) with the operator $\int_V (.) \Omega_k^*(\mathbf{x}) dv$ to yield the transformed algebraic system, in matrix form, as:

$$ (A + C) \{\bar{\psi}\} = \mu^2 E \{\bar{\psi}\} \qquad (13a) $$

with the elements of the coefficient matrices given by:

$$ a_{kl} = \frac{1}{M_l} \left[ \int_S \gamma_k (B - \hat{B}^*) \Omega_l(\mathbf{x}) ds + \int_V \Omega_k^*(\mathbf{x}) (L - \hat{L}^*) \Omega_l(\mathbf{x}) dv \right] \qquad (13b) $$

$$ c_{kl} = \lambda_k^2 \delta_{kl} \qquad (13c) $$

$$ e_{kl} = \frac{1}{M_l} \int_V w(\mathbf{x}) \Omega_k^*(\mathbf{x}) \Omega_l(\mathbf{x}) dv \qquad (13d) $$

$$ \mu = \{\mu_1, \mu_2, \dots\} \qquad (13e) $$

$$ \gamma_k = \frac{\Omega_k^*(\mathbf{x}) - n_l \hat{k}_{jl}(\mathbf{x}) \frac{\partial \tilde{\Omega}_k^*(\mathbf{x})}{\partial x_j}}{\hat{\alpha}(\mathbf{x}) + \hat{\beta}(\mathbf{x})} \qquad (13f) $$

Similarly, we proceed with the integral transformation of eq.(3a) with the operator $\int_V (.) \Omega_k(\mathbf{x}) dv$ to yield the transformed algebraic system for the eigenfunction $\psi_n^*(\mathbf{x})$, in matrix form, as:
Then, the transform-inverse pair is proposed for the desired eigenfunction, $\psi_n(\mathbf{x})$

as:

$$ \bar{\psi}_{n,k} = \int_V \hat{w}(\mathbf{x}) \Omega_k^*(\mathbf{x}) \psi_n(\mathbf{x}) dv \quad \text{transform} \quad (11a) $$

$$ \psi_n(\mathbf{x}) = \sum_{k=1}^{\infty} \frac{1}{M_k} \Omega_k(\mathbf{x}) \bar{\psi}_{n,k} \quad \text{inverse} \quad (11b) $$

Similarly, the transform-inverse pair for the adjoint eigenfunction, $\psi_n^*(\mathbf{x})$ becomes:

$$ \bar{\psi}_{n,k}^* = \int_V \hat{w}(\mathbf{x}) \Omega_k(\mathbf{x}) \psi_n^*(\mathbf{x}) dv \quad \text{transform} \quad (12a) $$

$$ \psi_n^*(\mathbf{x}) = \sum_{k=1}^{\infty} \frac{1}{M_k} \Omega_k^*(\mathbf{x}) \bar{\psi}_{n,k}^* \quad \text{inverse} \quad (12b) $$

Now, we proceed with the integral transformation of eq.(2a) with the operator $\int_V (.) \Omega_k^*(\mathbf{x}) dv$ to yield the transformed algebraic system, in matrix form, as:

$$ (A + C) \{\bar{\psi}\} = \mu^2 E \{\bar{\psi}\} \quad (13a) $$

with the elements of the coefficient matrices given by:

$$ a_{kl} = \frac{1}{M_l} \left[ \int_S \gamma_k (B - \hat{B}^*) \Omega_l(\mathbf{x}) ds + \int_V \Omega_k^*(\mathbf{x}) (L - \hat{L}^*) \Omega_l(\mathbf{x}) dv \right] \quad (13b) $$

$$ c_{kl} = \lambda_k^2 \delta_{kl} \quad (13c) $$

$$ e_{kl} = \frac{1}{M_l} \int_V w(\mathbf{x}) \Omega_k^*(\mathbf{x}) \Omega_l(\mathbf{x}) dv \quad (13d) $$

$$ \mu = \{\mu_1, \mu_2, \dots\} \quad (13e) $$

$$ \gamma_k = \frac{\Omega_k^*(\mathbf{x}) - n_k \hat{k}_{ji}(\mathbf{x}) \frac{\partial \tilde{\Omega}_k^*(\mathbf{x})}{\partial x_j}}{\hat{\alpha}(\mathbf{x}) + \hat{\beta}(\mathbf{x})} \quad (13f) $$

Similarly, we proceed with the integral transformation of eq.(3a) with the operator $\int_V (.) \Omega_k(\mathbf{x}) dv$ to yield the transformed algebraic system for the eigenfunction $\psi_n^*(\mathbf{x})$, in matrix form, as:
$$ (\mathbf{A}^* + \mathbf{C}) \{\bar{\psi}^*\} = \mu^2 \mathbf{E}^* \{\bar{\psi}^*\} \quad (14a) $$

with the elements of the coefficient matrices given by:

$$ a_{kl}^* = \frac{1}{M_l} \left[ \int_S \gamma_k^*(B^* - \hat{B})\Omega_l^*(\mathbf{x})ds + \int_V \Omega_k(\mathbf{x})(L^* - \hat{L})\Omega_l^*(\mathbf{x})dv \right] \quad (14b) $$

$$ e_{kl}^* = \frac{1}{M_l} \int_V w(\mathbf{x}) \Omega_l^*(\mathbf{x}) \Omega_k(\mathbf{x}) dv \quad (14c) $$

$$ \gamma_k^* = \frac{\Omega_k(\mathbf{x}) - n_i \hat{k}_{ij}(\mathbf{x}) \frac{\partial \tilde{\Omega}_k(\mathbf{x})}{\partial x_j}}{\hat{\alpha}(\mathbf{x}) + \hat{\beta}(\mathbf{x})} \quad (14d) $$

It is clear that if a self-adjoint auxiliary eigenvalue problem ($\hat{L}^* \equiv \hat{L}$, $\hat{B}^* \equiv \hat{B}$) is chosen to solve problems (2) and (3) above, there will be just one set of auxiliary eigenfunctions, $\Omega_k(\mathbf{x}) \equiv \Omega_k^*(\mathbf{x})$ and the matrix systems coefficients shall be simplified accordingly.

## GITT solution with self-adjoint eigenvalue problem

In the more traditional formalism of the GITT approach, one would usually consider the adoption of a classical self-adjoint Sturm-Liouville problem as the eigenfunction expansion base, while the cross-derivatives in the equation and the derivatives of the other coordinates in each boundary condition would be merged with the corresponding source terms. Thus, problem (1) would be rewritten as:

$$ w(\mathbf{x}) \frac{\partial T(\mathbf{x}, t)}{\partial t} + \tilde{L} T(\mathbf{x}, t) = \tilde{P}(\mathbf{x}, t, T), \quad \text{in } \mathbf{x} \in V, t > 0 \quad (15a) $$

with initial and boundary conditions given as

$$ T(\mathbf{x}, 0) = f(\mathbf{x}), \quad \mathbf{x} \in V \quad (15b) $$

$$ \tilde{B}T(\mathbf{x}, t) = \tilde{\phi}(\mathbf{x}, t, T), \quad \mathbf{x} \in S, t > 0 \quad (15c) $$

where the new equation and boundary operators corresponding to an orthotropic medium are written in tensor notation as

$$ \begin{aligned} \tilde{L} &= -\frac{\partial}{\partial x_i} \left( k_{ii}(\mathbf{x}) \frac{\partial}{\partial x_i} \right) + d(\mathbf{x}) \\ \tilde{B} &= \alpha(\mathbf{x}) + \beta(\mathbf{x}) n_i k_{ii}(\mathbf{x}) \frac{\partial}{\partial x_i} \end{aligned} \quad (15d,e) $$

and the source terms are redefined as
$$ (\mathbf{A}^* + \mathbf{C}) \{\bar{\psi}^*\} = \mu^2 \mathbf{E}^* \{\bar{\psi}^*\} \quad (14a) $$

with the elements of the coefficient matrices given by:

$$ a_{kl}^* = \frac{1}{M_l} \left[ \int_S \gamma_k^*(B^* - \hat{B})\Omega_l^*(\mathbf{x})ds + \int_V \Omega_k(\mathbf{x})(L^* - \hat{L})\Omega_l^*(\mathbf{x})dv \right] \quad (14b) $$

$$ e_{kl}^* = \frac{1}{M_l} \int_V w(\mathbf{x}) \Omega_l^*(\mathbf{x}) \Omega_k(\mathbf{x}) dv \quad (14c) $$

$$ \gamma_k^* = \frac{\Omega_k(\mathbf{x}) - n_i \hat{k}_{ij}(\mathbf{x}) \frac{\partial \tilde{\Omega}_k(\mathbf{x})}{\partial x_j}}{\hat{\alpha}(\mathbf{x}) + \hat{\beta}(\mathbf{x})} \quad (14d) $$

It is clear that if a self-adjoint auxiliary eigenvalue problem ($\hat{L}^* \equiv \hat{L}$, $\hat{B}^* \equiv \hat{B}$) is chosen to solve problems (2) and (3) above, there will be just one set of auxiliary eigenfunctions, $\Omega_k(\mathbf{x}) = \Omega_k^*(\mathbf{x})$ and the matrix systems coefficients shall be simplified accordingly.

## GITT solution with self-adjoint eigenvalue problem

In the more traditional formalism of the GITT approach, one would usually consider the adoption of a classical self-adjoint Sturm-Liouville problem as the eigenfunction expansion base, while the cross-derivatives in the equation and the derivatives of the other coordinates in each boundary condition would be merged with the corresponding source terms. Thus, problem (1) would be rewritten as:

$$ w(\mathbf{x}) \frac{\partial T(\mathbf{x}, t)}{\partial t} + \tilde{L} T(\mathbf{x}, t) = \tilde{P}(\mathbf{x}, t, T), \quad \text{in } \mathbf{x} \in V, t > 0 \quad (15a) $$

with initial and boundary conditions given as

$$ T(\mathbf{x}, 0) = f(\mathbf{x}), \quad \mathbf{x} \in V \quad (15b) $$

$$ \tilde{B}T(\mathbf{x}, t) = \tilde{\phi}(\mathbf{x}, t, T), \quad \mathbf{x} \in S, t > 0 \quad (15c) $$

where the new equation and boundary operators corresponding to an orthotropic medium are written in tensor notation as

$$ \begin{gathered} \tilde{L} = -\frac{\partial}{\partial x_i} \left( k_{ii}(\mathbf{x}) \frac{\partial}{\partial x_i} \right) + d(\mathbf{x}) \quad (15d,e) \\ \tilde{B} = \alpha(\mathbf{x}) + \beta(\mathbf{x}) n_i k_{ii}(\mathbf{x}) \frac{\partial}{\partial x_i} \end{gathered} $$

and the source terms are redefined as
$$ \tilde{P}(\mathbf{x}, t, T) = P(\mathbf{x}, t, T) + \frac{\partial}{\partial x_i} \left( (1 - \delta_{ij}) k_{ij}(\mathbf{x}) \frac{\partial}{\partial x_j} \right) \qquad (15f) $$

$$ \tilde{\phi}(\mathbf{x}, t, T) = \phi(\mathbf{x}, t, T) - \beta(\mathbf{x}) n_i k_{ij}(\mathbf{x}) (1 - \delta_{ij}) \frac{\partial T(\mathbf{x}, t)}{\partial x_j} \qquad (15g) $$

The self-adjoint eigenvalue problem is then chosen as:

$$ \tilde{L}\psi(\mathbf{x}) = \mu^2 w(\mathbf{x})\psi(\mathbf{x}), \quad \text{in } \mathbf{x} \in V \qquad (16a) $$

with boundary conditions given as

$$ \tilde{B}\psi(\mathbf{x}) = 0, \quad \mathbf{x} \in S \qquad (16b) $$

Following the formalism in the GITT approach [14,16], the integral transform-inverse pair is given by:

$$ \bar{T}_n(t) = \int_V w(\mathbf{x}) T(\mathbf{x}, t) \psi_n(\mathbf{x}) dv \quad \text{transform} \qquad (17a) $$

$$ T(\mathbf{x}, t) = \sum_{n=1}^{\infty} \frac{1}{N_n} \psi_n(\mathbf{x}) \bar{T}_n(t) \quad \text{inverse} \qquad (17b) $$

Then, the integral transformation of eq.(15a) is performed employing the operator $\int_V (.) \psi_n(\mathbf{x}) dv$, following the formalism in [14,16], to yield the transformed system below:

$$ \frac{d\bar{T}_n(t)}{dt} + \mu_n^2 \bar{T}_n(t) = \bar{g}_n(t, \bar{T}(t)), \quad n = 1, 2, \dots, t > 0 \qquad (18a) $$

$$ \bar{T}_n(0) = \bar{f}_n \qquad (18b) $$

where the transformed source terms and initial conditions are given by

$$ \bar{f}_n = \int_V w(\mathbf{x}) f(\mathbf{x}) \psi_n(\mathbf{x}) dv \qquad (18c) $$

$$ \bar{g}_n(t, \bar{T}(t)) = \int_S \tilde{\phi}(\mathbf{x}, t, T) \left[ \frac{\psi_n(\mathbf{x}) - n_i k_{ii}(\mathbf{x}) \frac{\partial \psi_n(\mathbf{x})}{\partial x_i}}{\alpha(\mathbf{x}) + \beta(\mathbf{x})} \right] ds + \int_V \tilde{P}(\mathbf{x}, t, T) \psi_n(\mathbf{x}) dv \qquad (18d) $$

Although fairly straightforward, the traditional GITT approach is here reviewed only for the sake of completeness, since it should lead to a slower convergence behavior in comparison to the here proposed extended methodology employing the adjoint eigenvalue problems that fully account for the general anisotropy. As can be seen from eq.(18d) above, the source terms in both the equation and boundary conditions carry all the information on the off-diagonal elements in the thermal conductivity tensor, and then compose the transformed source term, $\bar{g}_n(t, \bar{T}(t))$, responsible for the coupling of the
$$ \tilde{P}(\mathbf{x}, t, T) = P(\mathbf{x}, t, T) + \frac{\partial}{\partial x_i} \left( (1 - \delta_{ij}) k_{ij}(\mathbf{x}) \frac{\partial}{\partial x_j} \right) \qquad (15f) $$

$$ \tilde{\phi}(\mathbf{x}, t, T) = \phi(\mathbf{x}, t, T) - \beta(\mathbf{x}) n_i k_{ij}(\mathbf{x}) (1 - \delta_{ij}) \frac{\partial T(\mathbf{x}, t)}{\partial x_j} \qquad (15g) $$

The self-adjoint eigenvalue problem is then chosen as:

$$ \tilde{L}\psi(\mathbf{x}) = \mu^2 w(\mathbf{x})\psi(\mathbf{x}), \quad \text{in } \mathbf{x} \in V \qquad (16a) $$

with boundary conditions given as

$$ \tilde{B}\psi(\mathbf{x}) = 0, \quad \mathbf{x} \in S \qquad (16b) $$

Following the formalism in the GITT approach [14,16], the integral transform-inverse pair is given by:

$$ \bar{T}_n(t) = \int_V w(\mathbf{x}) T(\mathbf{x}, t) \psi_n(\mathbf{x}) dv \quad \text{transform} \qquad (17a) $$

$$ T(\mathbf{x}, t) = \sum_{n=1}^{\infty} \frac{1}{N_n} \psi_n(\mathbf{x}) \bar{T}_n(t) \quad \text{inverse} \qquad (17b) $$

Then, the integral transformation of eq.(15a) is performed employing the operator $\int_V (.) \psi_n(\mathbf{x}) dv$, following the formalism in [14,16], to yield the transformed system below:

$$ \frac{d\bar{T}_n(t)}{dt} + \mu_n^2 \bar{T}_n(t) = \bar{g}_n(t, \bar{T}(t)), \quad n = 1, 2, \dots, t > 0 \qquad (18a) $$

$$ \bar{T}_n(0) = \bar{f}_n \qquad (18b) $$

where the transformed source terms and initial conditions are given by

$$ \bar{f}_n = \int_V w(\mathbf{x}) f(\mathbf{x}) \psi_n(\mathbf{x}) dv \qquad (18c) $$

$$ \bar{g}_n(t, \bar{T}(t)) = \int_S \tilde{\phi}(\mathbf{x}, t, T) \left[ \frac{\psi_n(\mathbf{x}) - n_i k_{ii}(\mathbf{x}) \frac{\partial \psi_n(\mathbf{x})}{\partial x_i}}{\alpha(\mathbf{x}) + \beta(\mathbf{x})} \right] ds + \int_V \tilde{P}(\mathbf{x}, t, T) \psi_n(\mathbf{x}) dv \qquad (18d) $$

Although fairly straightforward, the traditional GITT approach is here reviewed only for the sake of completeness, since it should lead to a slower convergence behavior in comparison to the here proposed extended methodology employing the adjoint eigenvalue problems that fully account for the general anisotropy. As can be seen from eq.(18d) above, the source terms in both the equation and boundary conditions carry all the information on the off-diagonal elements in the thermal conductivity tensor, and then compose the transformed source term, $\bar{g}_n(t, \bar{T}(t))$, responsible for the coupling of the
transformed potentials and the slower convergence behavior. Analytical filtering and/or integral balance strategies [18,21,23] can in any case be employed to reduce the importance of these modified source terms. Also, the eigenvalue problem (16) can be readily solved for any arbitrary space variation of the coefficients following the same GITT procedure for eigenvalue problems above discussed [25-27], based on a simpler auxiliary self-adjoint eigenvalue problem.

## TEST CASE

To offer further details in the solution procedure, a two-dimensional transient heat conduction problem is considered, represented by an anisotropic material inclusion in a rectangular isotropic substrate. The formulation already considers homogeneous boundary conditions, thus skipping the step on the proposition of a filtering solution for homogenization. The test problem is then formulated as:

$$ w(x, y) \frac{\partial T(x, y, t)}{\partial t} = \frac{\partial}{\partial x} \left[ k_{11} \frac{\partial T}{\partial x} + k_{12} \frac{\partial T}{\partial y} \right] + \frac{\partial}{\partial y} \left[ k_{21} \frac{\partial T}{\partial x} + k_{22} \frac{\partial T}{\partial y} \right] - d(x, y)T + P(x, y, t), \quad \text{in } x \in [x_0, x_1], y \in [y_0, y_1], t > 0 \quad (19a) $$

with initial and boundary conditions given as

$$ T(x, y, 0) = f(x, y), \quad \text{in } x \in [x_0, x_1], y \in [y_0, y_1] \quad (19b) $$

$$ \alpha_{x,k} T + (-1)^{k+1} \beta_{x,k} \left[ k_{11} \frac{\partial T}{\partial x} + k_{12} \frac{\partial T}{\partial y} \right] = 0, \quad x = x_k, k=0,1, t > 0 \quad (19c) $$

$$ \alpha_{y,k} T + (-1)^{k+1} \beta_{y,k} \left[ k_{21} \frac{\partial T}{\partial x} + k_{22} \frac{\partial T}{\partial y} \right] = 0, \quad y = y_k, k=0,1, t > 0 \quad (19d) $$

Following the formalism described in the previous section, the biorthogonal set of eigenvalue problems is given by:

$$ \frac{\partial}{\partial x} \left[ k_{11} \frac{\partial \psi}{\partial x} + k_{12} \frac{\partial \psi}{\partial y} \right] + \frac{\partial}{\partial y} \left[ k_{21} \frac{\partial \psi}{\partial x} + k_{22} \frac{\partial \psi}{\partial y} \right] + \left[ \mu^2 w(x, y) - d(x, y) \right] \psi(x, y) = 0, \quad \text{in } x \in [x_0, x_1], y \in [y_0, y_1] \quad (20a) $$

with boundary conditions
transformed potentials and the slower convergence behavior. Analytical filtering and/or integral balance strategies [18,21,23] can in any case be employed to reduce the importance of these modified source terms. Also, the eigenvalue problem (16) can be readily solved for any arbitrary space variation of the coefficients following the same GITT procedure for eigenvalue problems above discussed [25-27], based on a simpler auxiliary self-adjoint eigenvalue problem.

## TEST CASE

To offer further details in the solution procedure, a two-dimensional transient heat conduction problem is considered, represented by an anisotropic material inclusion in a rectangular isotropic substrate. The formulation already considers homogeneous boundary conditions, thus skipping the step on the proposition of a filtering solution for homogenization. The test problem is then formulated as:

$$ w(x, y) \frac{\partial T(x, y, t)}{\partial t} = \frac{\partial}{\partial x} \left[ k_{11} \frac{\partial T}{\partial x} + k_{12} \frac{\partial T}{\partial y} \right] + \frac{\partial}{\partial y} \left[ k_{21} \frac{\partial T}{\partial x} + k_{22} \frac{\partial T}{\partial y} \right] - d(x, y)T + P(x, y, t), \quad \text{in } x \in [x_0, x_1], y \in [y_0, y_1], t > 0 \quad (19a) $$

with initial and boundary conditions given as

$$ T(x, y, 0) = f(x, y), \quad \text{in } x \in [x_0, x_1], y \in [y_0, y_1] \quad (19b) $$

$$ \alpha_{x,k} T + (-1)^{k+1} \beta_{x,k} \left[ k_{11} \frac{\partial T}{\partial x} + k_{12} \frac{\partial T}{\partial y} \right] = 0, \quad x = x_k, k=0,1, t > 0 \quad (19c) $$

$$ \alpha_{y,k} T + (-1)^{k+1} \beta_{y,k} \left[ k_{21} \frac{\partial T}{\partial x} + k_{22} \frac{\partial T}{\partial y} \right] = 0, \quad y = y_k, k=0,1, t > 0 \quad (19d) $$

Following the formalism described in the previous section, the biorthogonal set of eigenvalue problems is given by:

$$ \frac{\partial}{\partial x} \left[ k_{11} \frac{\partial \psi}{\partial x} + k_{12} \frac{\partial \psi}{\partial y} \right] + \frac{\partial}{\partial y} \left[ k_{21} \frac{\partial \psi}{\partial x} + k_{22} \frac{\partial \psi}{\partial y} \right] + \left[ \mu^2 w(x, y) - d(x, y) \right] \psi(x, y) = 0, \quad \text{in } x \in [x_0, x_1], y \in [y_0, y_1] \quad (20a) $$

with boundary conditions
$$ \alpha_{x,k} \psi + (-1)^{k+1} \beta_{x,k} \left[ k_{11} \frac{\partial \psi}{\partial x} + k_{12} \frac{\partial \psi}{\partial y} \right] = 0, \quad x = x_k, k=0,1 \quad (20b) $$

$$ \alpha_{y,k} \psi + (-1)^{k+1} \beta_{y,k} \left[ k_{21} \frac{\partial \psi}{\partial x} + k_{22} \frac{\partial \psi}{\partial y} \right] = 0, \quad y = y_k, k=0,1 \quad (20c) $$

and

$$ \frac{\partial}{\partial x} \left[ k_{11} \frac{\partial \psi^*}{\partial x} + k_{21} \frac{\partial \psi^*}{\partial y} \right] + \frac{\partial}{\partial y} \left[ k_{12} \frac{\partial \psi^*}{\partial x} + k_{22} \frac{\partial \psi^*}{\partial y} \right] + [\mu^2 w(x, y) - d(x, y)] \psi^*(x, y) = 0, \\ \text{in } x \in [x_0, x_1], y \in [y_0, y_1] \quad (21a) $$

with boundary conditions

$$ \alpha_{x,k} \psi^* + (-1)^{k+1} \beta_{x,k} \left[ k_{11} \frac{\partial \psi^*}{\partial x} + k_{21} \frac{\partial \psi^*}{\partial y} \right] = 0, \quad x = x_k, k=0,1 \quad (21b) $$

$$ \alpha_{y,k} \psi^* + (-1)^{k+1} \beta_{y,k} \left[ k_{12} \frac{\partial \psi^*}{\partial x} + k_{22} \frac{\partial \psi^*}{\partial y} \right] = 0, \quad y = y_k, k=0,1 \quad (21c) $$

The orthogonality property, Eq.(4a), is then given as:

$$ \int_{y_0}^{y_1} \int_{x_0}^{x_1} w(x, y) \psi_n(x, y) \psi_m^*(x, y) dx dy = \delta_{n,m} N_n \quad (22a) $$

with the normalization integral

$$ N_n = \int_{y_0}^{y_1} \int_{x_0}^{x_1} w(x, y) \psi_n(x, y) \psi_n^*(x, y) dx dy \quad (22b) $$

The integral transform pair then becomes:

$$ \bar{T}_n(t) = \int_{y_0}^{y_1} \int_{x_0}^{x_1} w(x, y) T(x, y, t) \psi_n^*(x, y) dx dy \quad \text{transform} \quad (23a) $$

$$ T(x, y, t) = \sum_{n=1}^{\infty} \frac{1}{N_n} \psi_n(x, y) \bar{T}_n(t) \quad \text{inverse} \quad (23b) $$

The transformed system results in:

$$ \frac{d\bar{T}_n(t)}{dt} + \mu_n^2 \bar{T}_n(t) = \bar{g}_n(t), \quad n=1,2,..., t>0 \quad (24a) $$

$$ \bar{T}_n(0) = \bar{f}_n \quad (24b) $$

where the transformed source terms and initial conditions are given by

$$ \bar{f}_n = \int_{y_0}^{y_1} \int_{x_0}^{x_1} w(x, y) f(x, y) \psi_n^*(x, y) dx dy \quad (24c) $$

$$ \bar{g}_n(t) = \int_{y_0}^{y_1} \int_{x_0}^{x_1} P(x, y, t) \psi_n^*(x, y) dx dy \quad (24d) $$
$$ \alpha_{x,k} \psi + (-1)^{k+1} \beta_{x,k} \left[ k_{11} \frac{\partial \psi}{\partial x} + k_{12} \frac{\partial \psi}{\partial y} \right] = 0, \quad x = x_k, k=0,1 \quad (20b) $$

$$ \alpha_{y,k} \psi + (-1)^{k+1} \beta_{y,k} \left[ k_{21} \frac{\partial \psi}{\partial x} + k_{22} \frac{\partial \psi}{\partial y} \right] = 0, \quad y = y_k, k=0,1 \quad (20c) $$

and

$$ \frac{\partial}{\partial x} \left[ k_{11} \frac{\partial \psi^*}{\partial x} + k_{21} \frac{\partial \psi^*}{\partial y} \right] + \frac{\partial}{\partial y} \left[ k_{12} \frac{\partial \psi^*}{\partial x} + k_{22} \frac{\partial \psi^*}{\partial y} \right] + [\mu^2 w(x, y) - d(x, y)] \psi^*(x, y) = 0, \quad \text{in } x \in [x_0, x_1], y \in [y_0, y_1] \quad (21a) $$

with boundary conditions

$$ \alpha_{x,k} \psi^* + (-1)^{k+1} \beta_{x,k} \left[ k_{11} \frac{\partial \psi^*}{\partial x} + k_{21} \frac{\partial \psi^*}{\partial y} \right] = 0, \quad x = x_k, k=0,1 \quad (21b) $$

$$ \alpha_{y,k} \psi^* + (-1)^{k+1} \beta_{y,k} \left[ k_{12} \frac{\partial \psi^*}{\partial x} + k_{22} \frac{\partial \psi^*}{\partial y} \right] = 0, \quad y = y_k, k=0,1 \quad (21c) $$

The orthogonality property, Eq.(4a), is then given as:

$$ \int_{y_0}^{y_1} \int_{x_0}^{x_1} w(x, y) \psi_n(x, y) \psi_m^*(x, y) dx dy = \delta_{n,m} N_n \quad (22a) $$

with the normalization integral

$$ N_n = \int_{y_0}^{y_1} \int_{x_0}^{x_1} w(x, y) \psi_n(x, y) \psi_n^*(x, y) dx dy \quad (22b) $$

The integral transform pair then becomes:

$$ \bar{T}_n(t) = \int_{y_0}^{y_1} \int_{x_0}^{x_1} w(x, y) T(x, y, t) \psi_n^*(x, y) dx dy \quad \text{transform} \quad (23a) $$

$$ T(x, y, t) = \sum_{n=1}^{\infty} \frac{1}{N_n} \psi_n(x, y) \bar{T}_n(t) \quad \text{inverse} \quad (23b) $$

The transformed system results in:

$$ \frac{d\bar{T}_n(t)}{dt} + \mu_n^2 \bar{T}_n(t) = \bar{g}_n(t), \quad n=1,2,..., t>0 \quad (24a) $$

$$ \bar{T}_n(0) = \bar{f}_n \quad (24b) $$

where the transformed source terms and initial conditions are given by

$$ \bar{f}_n = \int_{y_0}^{y_1} \int_{x_0}^{x_1} w(x, y) f(x, y) \psi_n^*(x, y) dx dy \quad (24c) $$

$$ \bar{g}_n(t) = \int_{y_0}^{y_1} \int_{x_0}^{x_1} P(x, y, t) \psi_n^*(x, y) dx dy \quad (24d) $$
Equations (24a,b) are then readily solvable analytically for this linear problem, as given by Eq.(7).

With respect to the solution of the two eigenvalue problems, Eqs.(20a-c) and (21a-c), one fairly straightforward choice of auxiliary eigenvalue problems corresponds to the orthotropic medium with constant coefficients, when $\Omega_k^*(x, y) = \Omega_k(x, y)$, given as:

$$ \hat{k}_{11} \frac{\partial^2 \Omega}{\partial x^2} + \hat{k}_{22} \frac{\partial^2 \Omega}{\partial y^2} + [\lambda^2 \hat{w} - \hat{d}] \Omega(x, y) = 0, \quad \text{in } x \in [x_0, x_1], y \in [y_0, y_1] \quad (25a) $$

with boundary conditions

$$ \alpha_{x,k} \Omega + (-1)^{k+1} \beta_{x,k} \left[ \hat{k}_{11} \frac{\partial \Omega}{\partial x} \right] = 0, \quad x = x_k, k=0,1 \quad (25b) $$

$$ \alpha_{y,k} \Omega + (-1)^{k+1} \beta_{y,k} \left[ \hat{k}_{22} \frac{\partial \Omega}{\partial y} \right] = 0, \quad y = y_k, k=0,1 \quad (25c) $$

where $\hat{k}_{11}$, $\hat{k}_{22}$, $\hat{w}$, and $\hat{d}$ are characteristic constant values of the original coefficients. Upon independent variables transformation, the auxiliary problem can be readily rewritten as:

$$ \frac{\partial^2 \Omega}{\partial X^2} + \frac{\partial^2 \Omega}{\partial Y^2} + \gamma^2 \Omega(X, Y) = 0, \quad \text{in } X \in [X_0, X_1], Y \in [Y_0, Y_1] \quad (26a) $$

with boundary conditions

$$ \alpha_{x,k} \Omega + (-1)^{k+1} \hat{\beta}_{x,k} \frac{\partial \Omega}{\partial X} = 0, \quad X = X_k, k=0,1 \quad (26b) $$

$$ \alpha_{y,k} \Omega + (-1)^{k+1} \hat{\beta}_{y,k} \frac{\partial \Omega}{\partial Y} = 0, \quad Y = Y_k, k=0,1 \quad (26c) $$

where,

$$ X = \frac{x}{\hat{k}_{11}^{1/2}}; \quad Y = \frac{y}{\hat{k}_{22}^{1/2}}; \quad X_k = \frac{x_k}{\hat{k}_{11}^{1/2}}; \quad Y_k = \frac{y_k}{\hat{k}_{22}^{1/2}}; \gamma^2 = \lambda^2 \hat{w} - \hat{d}; \quad \hat{\beta}_{x,k} = \frac{\beta_{x,k}}{\hat{k}_{11}^{1/2}}; \quad \hat{\beta}_{y,k} = \frac{\beta_{y,k}}{\hat{k}_{22}^{1/2}} \quad (26d-j) $$

and separation of variables is directly employed to find the analytical solution of the auxiliary eigenfunctions and eigenvalues.

This auxiliary eigenvalue problem is then employed in solving for the biorthogonal set of eigenfunctions, as illustrated below for $\Psi_n(x, y)$:

$$ \bar{\psi}_{n,k} = \int_{y_0}^{y_1} \int_{x_0}^{x_1} \hat{w} \tilde{\Omega}_k(x, y) \psi_n(x, y) dx dy \quad \text{transform} \quad (27a) $$
Equations (24a,b) are then readily solvable analytically for this linear problem, as given by Eq.(7).

With respect to the solution of the two eigenvalue problems, Eqs.(20a-c) and (21a-c), one fairly straightforward choice of auxiliary eigenvalue problems corresponds to the orthotropic medium with constant coefficients, when $\Omega_k^*(x, y) = \Omega_k(x, y)$, given as:

$$ \hat{k}_{11} \frac{\partial^2 \Omega}{\partial x^2} + \hat{k}_{22} \frac{\partial^2 \Omega}{\partial y^2} + [\lambda^2 \hat{w} - \hat{d}] \Omega(x, y) = 0, \quad \text{in } x \in [x_0, x_1], y \in [y_0, y_1] \quad (25a) $$

with boundary conditions

$$ \alpha_{x,k} \Omega + (-1)^{k+1} \beta_{x,k} \left[ \hat{k}_{11} \frac{\partial \Omega}{\partial x} \right] = 0, \quad x = x_k, k=0,1 \quad (25b) $$

$$ \alpha_{y,k} \Omega + (-1)^{k+1} \beta_{y,k} \left[ \hat{k}_{22} \frac{\partial \Omega}{\partial y} \right] = 0, \quad y = y_k, k=0,1 \quad (25c) $$

where $\hat{k}_{11}$, $\hat{k}_{22}$, $\hat{w}$, and $\hat{d}$ are characteristic constant values of the original coefficients. Upon independent variables transformation, the auxiliary problem can be readily rewritten as:

$$ \frac{\partial^2 \Omega}{\partial X^2} + \frac{\partial^2 \Omega}{\partial Y^2} + \gamma^2 \Omega(X, Y) = 0, \quad \text{in } X \in [X_0, X_1], Y \in [Y_0, Y_1] \quad (26a) $$

with boundary conditions

$$ \alpha_{x,k} \Omega + (-1)^{k+1} \hat{\beta}_{x,k} \frac{\partial \Omega}{\partial X} = 0, \quad X = X_k, k=0,1 \quad (26b) $$

$$ \alpha_{y,k} \Omega + (-1)^{k+1} \hat{\beta}_{y,k} \frac{\partial \Omega}{\partial Y} = 0, \quad Y = Y_k, k=0,1 \quad (26c) $$

where,

$$ X = \frac{x}{\hat{k}_{11}^{1/2}}; \quad Y = \frac{y}{\hat{k}_{22}^{1/2}}; \quad X_k = \frac{x_k}{\hat{k}_{11}^{1/2}}; \quad Y_k = \frac{y_k}{\hat{k}_{22}^{1/2}}; \gamma^2 = \lambda^2 \hat{w} - \hat{d}; \quad \hat{\beta}_{x,k} = \frac{\beta_{x,k}}{\hat{k}_{11}^{1/2}}; \quad \hat{\beta}_{y,k} = \frac{\beta_{y,k}}{\hat{k}_{22}^{1/2}} \quad (26d-j) $$

and separation of variables is directly employed to find the analytical solution of the auxiliary eigenfunctions and eigenvalues.

This auxiliary eigenvalue problem is then employed in solving for the biorthogonal set of eigenfunctions, as illustrated below for $\Psi_n(x, y)$:

$$ \bar{\psi}_{n,k} = \int_{y_0}^{y_1} \int_{x_0}^{x_1} \hat{w} \tilde{\Omega}_k(x, y) \psi_n(x, y) dx dy \quad \text{transform} \quad (27a) $$
$$ \psi_n(x, y) = \sum_{k=1}^{\infty} \tilde{\Omega}_k(x, y) \bar{\psi}_{n,k} \quad \text{inverse} \qquad (27b) $$

Also, it has been preferred to express the integral transform pair using the normalized auxiliary eigenfunction, given as

$$ \tilde{\Omega}_k(x, y) = \frac{1}{M_k^{1/2}} \Omega_k(x, y) \qquad (27c) $$

with normalization integral

$$ M_k = \int_{y_0}^{y_1} \int_{x_0}^{x_1} \hat{w} \Omega_k^2(x, y) dx dy \qquad (27d) $$

Following the integral transformation procedure as applied first to problem (20), for instance, the resulting matrix eigenvalue problem is written in matrix form as in eq.(13a), with the elements of the coefficient matrices given by:

$$ a_{nm} = a_{1,nm} - (a_{2,nm} + a_{3,nm}) + a_{4,nm} - (a_{5,nm} + a_{6,nm}) \qquad (28a) $$

$$ c_{nm} = - \int_{y_0}^{y_1} \int_{x_0}^{x_1} d(x, y) \tilde{\Omega}_m(x, y) \tilde{\Omega}_n(x, y) dx dy \qquad (28b) $$

$$ e_{nm} = \int_{y_0}^{y_1} \int_{x_0}^{x_1} w(x, y) \tilde{\Omega}_m(x, y) \tilde{\Omega}_n(x, y) dx dy \qquad (28c) $$

where,

$$ a_{1,nm} = \int_{y_0}^{y_1} \left[ \left( k_{11} \frac{\partial \tilde{\Omega}_m}{\partial x} + k_{12} \frac{\partial \tilde{\Omega}_m}{\partial y} \right) \tilde{\Omega}_n \Bigg|_{x=x_1} - \left( k_{11} \frac{\partial \tilde{\Omega}_m}{\partial x} + k_{12} \frac{\partial \tilde{\Omega}_m}{\partial y} \right) \tilde{\Omega}_n \Bigg|_{x=x_0} \right] dy \qquad (28d) $$

$$ a_{2,nm} = \int_{y_0}^{y_1} \int_{x_0}^{x_1} k_{11} \frac{\partial \tilde{\Omega}_m}{\partial x} \frac{\partial \tilde{\Omega}_n}{\partial x} dx dy \qquad (28e) $$

$$ a_{3,nm} = \int_{y_0}^{y_1} \int_{x_0}^{x_1} k_{12} \frac{\partial \tilde{\Omega}_m}{\partial y} \frac{\partial \tilde{\Omega}_n}{\partial x} dx dy \qquad (27f) $$

$$ a_{4,nm} = \int_{x_0}^{x_1} \left[ \left( k_{21} \frac{\partial \tilde{\Omega}_m}{\partial x} + k_{12} \frac{\partial \tilde{\Omega}_m}{\partial y} \right) \tilde{\Omega}_n \Bigg|_{y=y_1} - \left( k_{21} \frac{\partial \tilde{\Omega}_m}{\partial x} + k_{22} \frac{\partial \tilde{\Omega}_m}{\partial y} \right) \tilde{\Omega}_n \Bigg|_{y=y_0} \right] dx \qquad (27g) $$

$$ a_{5,nm} = \int_{y_0}^{y_1} \int_{x_0}^{x_1} k_{21} \frac{\partial \tilde{\Omega}_m}{\partial x} \frac{\partial \tilde{\Omega}_n}{\partial y} dx dy \qquad (27h) $$

$$ a_{6,nm} = \int_{y_0}^{y_1} \int_{x_0}^{x_1} k_{22} \frac{\partial \tilde{\Omega}_m}{\partial y} \frac{\partial \tilde{\Omega}_n}{\partial y} dx dy \qquad (27i) $$
$$ \psi_n(x, y) = \sum_{k=1}^{\infty} \tilde{\Omega}_k(x, y) \bar{\psi}_{n,k} \quad \text{inverse} \qquad (27b) $$

Also, it has been preferred to express the integral transform pair using the normalized auxiliary eigenfunction, given as

$$ \tilde{\Omega}_k(x, y) = \frac{1}{M_k^{1/2}} \Omega_k(x, y) \qquad (27c) $$

with normalization integral

$$ M_k = \int_{y_0}^{y_1} \int_{x_0}^{x_1} \hat{w} \Omega_k^2(x, y) dx dy \qquad (27d) $$

Following the integral transformation procedure as applied first to problem (20), for instance, the resulting matrix eigenvalue problem is written in matrix form as in eq.(13a), with the elements of the coefficient matrices given by:

$$ a_{nm} = a_{1,nm} - (a_{2,nm} + a_{3,nm}) + a_{4,nm} - (a_{5,nm} + a_{6,nm}) \qquad (28a) $$

$$ c_{nm} = - \int_{y_0}^{y_1} \int_{x_0}^{x_1} d(x, y) \tilde{\Omega}_m(x, y) \tilde{\Omega}_n(x, y) dx dy \qquad (28b) $$

$$ e_{nm} = \int_{y_0}^{y_1} \int_{x_0}^{x_1} w(x, y) \tilde{\Omega}_m(x, y) \tilde{\Omega}_n(x, y) dx dy \qquad (28c) $$

where,

$$ a_{1,nm} = \int_{y_0}^{y_1} \left[ \left( k_{11} \frac{\partial \tilde{\Omega}_m}{\partial x} + k_{12} \frac{\partial \tilde{\Omega}_m}{\partial y} \right) \tilde{\Omega}_n \Bigg|_{x=x_1} - \left( k_{11} \frac{\partial \tilde{\Omega}_m}{\partial x} + k_{12} \frac{\partial \tilde{\Omega}_m}{\partial y} \right) \tilde{\Omega}_n \Bigg|_{x=x_0} \right] dy \qquad (28d) $$

$$ a_{2,nm} = \int_{y_0}^{y_1} \int_{x_0}^{x_1} k_{11} \frac{\partial \tilde{\Omega}_m}{\partial x} \frac{\partial \tilde{\Omega}_n}{\partial x} dx dy \qquad (28e) $$

$$ a_{3,nm} = \int_{y_0}^{y_1} \int_{x_0}^{x_1} k_{12} \frac{\partial \tilde{\Omega}_m}{\partial y} \frac{\partial \tilde{\Omega}_n}{\partial x} dx dy \qquad (27f) $$

$$ a_{4,nm} = \int_{x_0}^{x_1} \left[ \left( k_{21} \frac{\partial \tilde{\Omega}_m}{\partial x} + k_{12} \frac{\partial \tilde{\Omega}_m}{\partial y} \right) \tilde{\Omega}_n \Bigg|_{y=y_1} - \left( k_{21} \frac{\partial \tilde{\Omega}_m}{\partial x} + k_{22} \frac{\partial \tilde{\Omega}_m}{\partial y} \right) \tilde{\Omega}_n \Bigg|_{y=y_0} \right] dx \qquad (27g) $$

$$ a_{5,nm} = \int_{y_0}^{y_1} \int_{x_0}^{x_1} k_{21} \frac{\partial \tilde{\Omega}_m}{\partial x} \frac{\partial \tilde{\Omega}_n}{\partial y} dx dy \qquad (27h) $$

$$ a_{6,nm} = \int_{y_0}^{y_1} \int_{x_0}^{x_1} k_{22} \frac{\partial \tilde{\Omega}_m}{\partial y} \frac{\partial \tilde{\Omega}_n}{\partial y} dx dy \qquad (27i) $$
RESULTS AND DISCUSSION

A symbolic-numerical algorithm was constructed in the *Mathematica* v.13.2 software system, for the GITT solution of problem (19) with arbitrarily two-dimensional spatial variable coefficients. In order to report some numerical results, the required data was provided for the test case, represented by an isotropic substrate given by a square plate of side *a*=*b*=0.08 m, thermal conductivity *k*<sub>s</sub>=3.6 W/m°C, specific mass ρ<sub>s</sub>=3100 kg/m<sup>3</sup>, and specific heat *c*<sub>p,s</sub>=1100 J/kg°C. An anisotropic material inclusion is considered as a square insert in the middle of the substrate plate, with sides *l*<sub>x</sub>=*l*<sub>y</sub>=0.015m, thermal conductivity components *k*<sub>11</sub>= *k*<sub>22</sub>=26 W/m°C , *k*<sub>12</sub>= *k*<sub>21</sub>=8.667 W/m°C , specific mass ρ<sub>c</sub>=8666 kg/m<sup>3</sup>, and specific heat *c*<sub>p,c</sub>=343 J/kg°C. There is no heat generation source (*P*=0), no linear heat dissipation (*d*=0), and the initial condition is an uniform temperature field, *f*=*T*<sub>in</sub>=100°C. The boundary conditions at *x*=*y*=0 are insulated boundaries, α<sub>x,0</sub> = α<sub>y,0</sub> =0 and β<sub>x,0</sub>=β<sub>y,0</sub>=1, and at the boundaries *x*=*y*=*a*, prescribed temperatures (*T*<sub>w</sub>=0°C) or α<sub>x,1</sub> = α<sub>y,1</sub>=1 and β<sub>x,1</sub>=β<sub>y,1</sub>=0. The test case is simulated up to a final time value of *t*<sub>f</sub>=3600 seconds. The proposed heterogeneous test case then involves an abrupt variation of thermophysical properties at the interfaces of the substrate and central inclusion regions.

The first step in the solution procedure is the integral transform solution of the eigenvalue problems, eqs. (20,21), based on the analytical solution of the auxiliary problems, eqs. (26), taking the properties of the substrate as characteristic values. Considering different truncation orders in the inverse formula, eq.(27b), one may inspect the convergence behavior in the eigenfunction expansion proposal. Table I below provides numerical results for a few selected eigenvalues of problem (19), including the five most fundamental modes (μ₁ to μ₅), for truncation orders that vary from *M*=50 to 600. The last column provides a set of converged results to the four significant digits provided (plus or minus one in the last significant digit). One may observe that fairly low truncation orders already provide two or three converged digits, and that the estimate will degenerate for eigenvalue orders close to the truncation order, such as the result for μ₅₀ with *M*=50.
RESULTS AND DISCUSSION

A symbolic-numerical algorithm was constructed in the *Mathematica* v.13.2 software system, for the GITT solution of problem (19) with arbitrarily two-dimensional spatial variable coefficients. In order to report some numerical results, the required data was provided for the test case, represented by an isotropic substrate given by a square plate of side *a*=*b*=0.08 m, thermal conductivity *k*<sub>s</sub>=3.6 W/m°C, specific mass ρ<sub>s</sub>=3100 kg/m<sup>3</sup>, and specific heat *c*<sub>p,s</sub>=1100 J/kg°C. An anisotropic material inclusion is considered as a square insert in the middle of the substrate plate, with sides *l*<sub>x</sub>=*l*<sub>y</sub>=0.015m, thermal conductivity components *k*<sub>11</sub>= *k*<sub>22</sub>=26 W/m°C , *k*<sub>12</sub>= *k*<sub>21</sub>=8.667 W/m°C , specific mass ρ<sub>c</sub>=8666 kg/m<sup>3</sup>, and specific heat *c*<sub>p,c</sub>=343 J/kg°C. There is no heat generation source (*P*=0), no linear heat dissipation (*d*=0), and the initial condition is an uniform temperature field, *f*=*T*<sub>in</sub>=100°C. The boundary conditions at *x*=*y*=0 are insulated boundaries, α<sub>x,0</sub> = α<sub>y,0</sub> =0 and β<sub>x,0</sub>=β<sub>y,0</sub>=1, and at the boundaries *x*=*y*=*a*, prescribed temperatures (*T*<sub>w</sub>=0°C) or α<sub>x,1</sub> = α<sub>y,1</sub>=1 and β<sub>x,1</sub>=β<sub>y,1</sub>=0. The test case is simulated up to a final time value of *t*<sub>f</sub>=3600 seconds. The proposed heterogeneous test case then involves an abrupt variation of thermophysical properties at the interfaces of the substrate and central inclusion regions.

The first step in the solution procedure is the integral transform solution of the eigenvalue problems, eqs. (20,21), based on the analytical solution of the auxiliary problems, eqs. (26), taking the properties of the substrate as characteristic values. Considering different truncation orders in the inverse formula, eq.(27b), one may inspect the convergence behavior in the eigenfunction expansion proposal. Table I below provides numerical results for a few selected eigenvalues of problem (19), including the five most fundamental modes (μ₁ to μ₅), for truncation orders that vary from M=50 to 600. The last column provides a set of converged results to the four significant digits provided (plus or minus one in the last significant digit). One may observe that fairly low truncation orders already provide two or three converged digits, and that the estimate will degenerate for eigenvalue orders close to the truncation order, such as the result for μ₅₀ with M=50.
<table><caption>Table I – Convergence of selected eigenvalues in terms of truncation order M.</caption><thead><tr><th>M</th><th>50</th><th>100</th><th>200</th><th>300</th><th>400</th><th>500</th><th>600</th></tr></thead><tbody><tr><th>μ</th><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>1</th><td>0.02981</td><td>0.02975</td><td>0.02963</td><td>0.02961</td><td>0.02959</td><td>0.02957</td><td>0.02956</td></tr><tr><th>2</th><td>0.06521</td><td>0.06513</td><td>0.06502</td><td>0.06500</td><td>0.06498</td><td>0.06496</td><td>0.06495</td></tr><tr><th>3</th><td>0.06671</td><td>0.06660</td><td>0.06639</td><td>0.06637</td><td>0.06632</td><td>0.06629</td><td>0.06627</td></tr><tr><th>4</th><td>0.08878</td><td>0.08863</td><td>0.08837</td><td>0.08833</td><td>0.08828</td><td>0.08823</td><td>0.08821</td></tr><tr><th>5</th><td>0.1047</td><td>0.1045</td><td>0.1044</td><td>0.1044</td><td>0.1043</td><td>0.1043</td><td>0.1043</td></tr><tr><th>10</th><td>0.1470</td><td>0.1468</td><td>0.1466</td><td>0.1465</td><td>0.1465</td><td>0.1464</td><td>0.1464</td></tr><tr><th>20</th><td>0.2063</td><td>0.2057</td><td>0.2051</td><td>0.2050</td><td>0.2049</td><td>0.2048</td><td>0.2048</td></tr><tr><th>30</th><td>0.2595</td><td>0.2562</td><td>0.2551</td><td>0.2548</td><td>0.2546</td><td>0.2544</td><td>0.2543</td></tr><tr><th>40</th><td>0.3006</td><td>0.2924</td><td>0.2919</td><td>0.2919</td><td>0.2918</td><td>0.2917</td><td>0.2917</td></tr><tr><th>50</th><td>0.4976</td><td>0.3311</td><td>0.3301</td><td>0.3296</td><td>0.3292</td><td>0.3291</td><td>0.3289</td></tr></tbody></table>

Next, one may examine the convergence behavior of the temperature field proposed expansion, eq.(23b), in terms of the truncation order of the temperature eigenfunction expansion (N). Taking a fixed truncation order for the eigenvalue problem expansion (M=600), Table II provides temperature values at different positions within both the substrate and anisotropic insert and different time values, with increasing truncation orders in eq.(23b). Clearly, with fairly low truncation orders below N=50, the temperature values are already fully converged to the six significant digits provided, and even to the fifth digit with N<25, showing that the choice of a biorthogonal eigenvalue problems set indeed yields a rapidly converging temperature expansion.

<table><caption>Table II – Convergence of transient temperature in terms of temperature expansion truncation order (N) with M=600 terms for eigenvalue problem.</caption><thead><tr><th>N<br>(x, y, t)</th><th>10</th><th>25</th><th>50</th><th>75</th><th>100</th></tr></thead><tbody><tr><th>(0.02,0.02,180)</th><td>99.1125</td><td>99.3510</td><td>99.3513</td><td>99.3513</td><td>99.3513</td></tr><tr><th>(0.04,0.04,180)</th><td>90.0926</td><td>90.1142</td><td>90.1144</td><td>90.1144</td><td>90.1144</td></tr><tr><th>(0.06,0.06,180)</th><td>48.3814</td><td>48.5623</td><td>48.5616</td><td>48.5616</td><td>48.5616</td></tr><tr><th>(0.02,0.02,360)</th><td>92.6351</td><td>92.6405</td><td>92.6405</td><td>92.6405</td><td>92.6405</td></tr></tbody></table>
<table><caption>Table I – Convergence of selected eigenvalues in terms of truncation order M.</caption><thead><tr><th>M</th><th>50</th><th>100</th><th>200</th><th>300</th><th>400</th><th>500</th><th>600</th></tr></thead><tbody><tr><th>μ</th><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>1</th><td>0.02981</td><td>0.02975</td><td>0.02963</td><td>0.02961</td><td>0.02959</td><td>0.02957</td><td>0.02956</td></tr><tr><th>2</th><td>0.06521</td><td>0.06513</td><td>0.06502</td><td>0.06500</td><td>0.06498</td><td>0.06496</td><td>0.06495</td></tr><tr><th>3</th><td>0.06671</td><td>0.06660</td><td>0.06639</td><td>0.06637</td><td>0.06632</td><td>0.06629</td><td>0.06627</td></tr><tr><th>4</th><td>0.08878</td><td>0.08863</td><td>0.08837</td><td>0.08833</td><td>0.08828</td><td>0.08823</td><td>0.08821</td></tr><tr><th>5</th><td>0.1047</td><td>0.1045</td><td>0.1044</td><td>0.1044</td><td>0.1043</td><td>0.1043</td><td>0.1043</td></tr><tr><th>10</th><td>0.1470</td><td>0.1468</td><td>0.1466</td><td>0.1465</td><td>0.1465</td><td>0.1464</td><td>0.1464</td></tr><tr><th>20</th><td>0.2063</td><td>0.2057</td><td>0.2051</td><td>0.2050</td><td>0.2049</td><td>0.2048</td><td>0.2048</td></tr><tr><th>30</th><td>0.2595</td><td>0.2562</td><td>0.2551</td><td>0.2548</td><td>0.2546</td><td>0.2544</td><td>0.2543</td></tr><tr><th>40</th><td>0.3006</td><td>0.2924</td><td>0.2919</td><td>0.2919</td><td>0.2918</td><td>0.2917</td><td>0.2917</td></tr><tr><th>50</th><td>0.4976</td><td>0.3311</td><td>0.3301</td><td>0.3296</td><td>0.3292</td><td>0.3291</td><td>0.3289</td></tr></tbody></table>

Next, one may examine the convergence behavior of the temperature field proposed expansion, eq.(23b), in terms of the truncation order of the temperature eigenfunction expansion (N). Taking a fixed truncation order for the eigenvalue problem expansion (M=600), Table II provides temperature values at different positions within both the substrate and anisotropic insert and different time values, with increasing truncation orders in eq.(23b). Clearly, with fairly low truncation orders below N=50, the temperature values are already fully converged to the six significant digits provided, and even to the fifth digit with N<25, showing that the choice of a biorthogonal eigenvalue problems set indeed yields a rapidly converging temperature expansion.

<table><caption>Table II – Convergence of transient temperature in terms of temperature expansion truncation order (N) with M=600 terms for eigenvalue problem.</caption><thead><tr><th>N<br>(x, y, t)</th><th>10</th><th>25</th><th>50</th><th>75</th><th>100</th></tr></thead><tbody><tr><th>(0.02,0.02,180)</th><td>99.1125</td><td>99.3510</td><td>99.3513</td><td>99.3513</td><td>99.3513</td></tr><tr><th>(0.04,0.04,180)</th><td>90.0926</td><td>90.1142</td><td>90.1144</td><td>90.1144</td><td>90.1144</td></tr><tr><th>(0.06,0.06,180)</th><td>48.3814</td><td>48.5623</td><td>48.5616</td><td>48.5616</td><td>48.5616</td></tr><tr><th>(0.02,0.02,360)</th><td>92.6351</td><td>92.6405</td><td>92.6405</td><td>92.6405</td><td>92.6405</td></tr></tbody></table>
<table><tr><td>(0.04,0.04,360)</td><td>70.6628</td><td>70.6636</td><td>70.6636</td><td>70.6636</td><td>70.6636</td></tr><tr><td>(0.06,0.06,360)</td><td>29.5010</td><td>29.5041</td><td>29.5041</td><td>29.5041</td><td>29.5041</td></tr></table>

In order to offer an illustration of the GITT usefulness in benchmarking for purely numerical methodologies and codes, the test case was also solved employing the well-known COMSOL v.5.6 finite element code [34], considering different levels of mesh refinement as automatically provided by the software package. Table III below provides a summary of the different mesh sizes considered, including the adopted time step in each illustrative run. Figure 1 below provides a sample of the meshing for the Extra Fine case.

<table><caption>Table III – Different mesh sizes and time steps considered in COMSOL solutions.</caption><thead><tr><td>Name</td><td>Number of elements</td><td>Δt (s)</td></tr></thead><tbody><tr><td>Coarse</td><td>282</td><td>1</td></tr><tr><td>Normal</td><td>634</td><td>0.125</td></tr><tr><td>Fine</td><td>1004</td><td>0.084</td></tr><tr><td>Finer</td><td>1950</td><td>0.066</td></tr><tr><td>Extra Fine</td><td>6680</td><td>0.046</td></tr></tbody></table>

<figure><img src="image_2.png" /></figure>
<table><tr><td>(0.04,0.04,360)</td><td>70.6628</td><td>70.6636</td><td>70.6636</td><td>70.6636</td><td>70.6636</td></tr><tr><td>(0.06,0.06,360)</td><td>29.5010</td><td>29.5041</td><td>29.5041</td><td>29.5041</td><td>29.5041</td></tr></table>

In order to offer an illustration of the GITT usefulness in benchmarking for purely numerical methodologies and codes, the test case was also solved employing the well-known COMSOL v.5.6 finite element code [34], considering different levels of mesh refinement as automatically provided by the software package. Table III below provides a summary of the different mesh sizes considered, including the adopted time step in each illustrative run. Figure 1 below provides a sample of the meshing for the Extra Fine case.

<table><caption>Table III – Different mesh sizes and time steps considered in COMSOL solutions.</caption><thead><tr><td>Name</td><td>Number of elements</td><td>Δt (s)</td></tr></thead><tbody><tr><td>Coarse</td><td>282</td><td>1</td></tr><tr><td>Normal</td><td>634</td><td>0.125</td></tr><tr><td>Fine</td><td>1004</td><td>0.084</td></tr><tr><td>Finer</td><td>1950</td><td>0.066</td></tr><tr><td>Extra Fine</td><td>6680</td><td>0.046</td></tr></tbody></table>

<figure><img src="image_2.png" /></figure>
**Figure 1 – Finite element meshing for the Extra Fine mesh refinement option in COMSOL v.5.6 [34].**

The convergence of the eigenvalues and corresponding eigenfunctions, briefly illustrated in Table I, has a direct influence in the temperature field convergence as now analyzed in Table IV, where the truncation order *M* is increased from 100 to 600 terms, while keeping *N*=*M* in the temperature expansion. The results for *M*=600 are fully converged to plus or minus one in the fourth significant digit, and the corresponding column provides a set of benchmark results for reference purposes. Clearly, the convergence rates of the two expansions, for temperatures and eigenfunctions, have to be analyzed simultaneously, so as to surely yield fully converged results for the temperature distribution. Also shown are the results for the Extra Fine mesh refinement option of Table III from COMSOL [34], for the same positions and time values.

<table><caption>**Table IV – Convergence of transient temperature at different positions with *M* terms in the eigenvalue problem expansion truncation order (*N*=*M*).**</caption><thead><tr><th>M<br>(x, y, t)</th><th>100</th><th>200</th><th>300</th><th>400</th><th>500</th><th>600</th><th>COMSOL</th></tr></thead><tbody><tr><td>(0.02,0.02,180)</td><td>99.24</td><td>99.33</td><td>99.27</td><td>99.36</td><td>99.37</td><td>99.35</td><td>98.82</td></tr><tr><td>(0.04,0.04,180)</td><td>89.81</td><td>89.95</td><td>90.11</td><td>90.14</td><td>90.11</td><td>90.11</td><td>89.96</td></tr><tr><td>(0.06,0.06,180)</td><td>48.79</td><td>48.57</td><td>48.59</td><td>48.58</td><td>48.57</td><td>48.56</td><td>51.69</td></tr><tr><td>(0.02,0.02,360)</td><td>92.26</td><td>92.52</td><td>92.48</td><td>92.60</td><td>92.66</td><td>92.64</td><td>92.03</td></tr><tr><td>(0.04,0.04,360)</td><td>70.383</td><td>70.50</td><td>70.65</td><td>70.69</td><td>70.66</td><td>70.66</td><td>72.34</td></tr><tr><td>(0.06,0.06,360)</td><td>30.00</td><td>29.60</td><td>29.60</td><td>29.55</td><td>29.52</td><td>29.50</td><td>31.98</td></tr><tr><td>(0.02,0.02,720)</td><td>71.16</td><td>71.60</td><td>71.60</td><td>71.74</td><td>71.83</td><td>71.83</td><td>72.79</td></tr><tr><td>(0.04,0.04,720)</td><td>46.32</td><td>46.41</td><td>46.50</td><td>46.53</td><td>46.52</td><td>46.52</td><td>48.25</td></tr><tr><td>(0.06,0.06,720)</td><td>17.50</td><td>17.11</td><td>17.10</td><td>17.04</td><td>17.01</td><td>16.99</td><td>17.75</td></tr><tr><td>(0.02,0.02,1800)</td><td>27.98</td><td>28.38</td><td>28.41</td><td>28.52</td><td>28.59</td><td>28.61</td><td>29.06</td></tr><tr><td>(0.04,0.04,1800)</td><td>17.02</td><td>17.17</td><td>17.22</td><td>17.25</td><td>17.27</td><td>17.28</td><td>17.12</td></tr><tr><td>(0.06,0.06,1800)</td><td>6.160</td><td>6.036</td><td>6.036</td><td>6.016</td><td>6.012</td><td>6.005</td><td>5.76</td></tr></tbody></table>

The COMSOL numerical solution is also analyzed in terms of the dimensionless temperature, $\theta=(T_{in}-T)/T_{in}$, relative deviations with respect to the converged GITT results, Table V. The deviations are quite satisfactory, being throughout most of the time and space domains below 5%.
Figure 1 – Finite element meshing for the Extra Fine mesh refinement option in COMSOL v.5.6 [34].

The convergence of the eigenvalues and corresponding eigenfunctions, briefly illustrated in Table I, has a direct influence in the temperature field convergence as now analyzed in Table IV, where the truncation order *M* is increased from 100 to 600 terms, while keeping *N*=*M* in the temperature expansion. The results for *M*=600 are fully converged to plus or minus one in the fourth significant digit, and the corresponding column provides a set of benchmark results for reference purposes. Clearly, the convergence rates of the two expansions, for temperatures and eigenfunctions, have to be analyzed simultaneously, so as to surely yield fully converged results for the temperature distribution. Also shown are the results for the Extra Fine mesh refinement option of Table III from COMSOL [34], for the same positions and time values.

<table><caption>Table IV – Convergence of transient temperature at different positions with *M* terms in the eigenvalue problem expansion truncation order (*N*=*M*).</caption><thead><tr><th>M<br>(x, y, t)</th><th>100</th><th>200</th><th>300</th><th>400</th><th>500</th><th>600</th><th>COMSOL</th></tr></thead><tbody><tr><td>(0.02,0.02,180)</td><td>99.24</td><td>99.33</td><td>99.27</td><td>99.36</td><td>99.37</td><td>99.35</td><td>98.82</td></tr><tr><td>(0.04,0.04,180)</td><td>89.81</td><td>89.95</td><td>90.11</td><td>90.14</td><td>90.11</td><td>90.11</td><td>89.96</td></tr><tr><td>(0.06,0.06,180)</td><td>48.79</td><td>48.57</td><td>48.59</td><td>48.58</td><td>48.57</td><td>48.56</td><td>51.69</td></tr><tr><td>(0.02,0.02,360)</td><td>92.26</td><td>92.52</td><td>92.48</td><td>92.60</td><td>92.66</td><td>92.64</td><td>92.03</td></tr><tr><td>(0.04,0.04,360)</td><td>70.383</td><td>70.50</td><td>70.65</td><td>70.69</td><td>70.66</td><td>70.66</td><td>72.34</td></tr><tr><td>(0.06,0.06,360)</td><td>30.00</td><td>29.60</td><td>29.60</td><td>29.55</td><td>29.52</td><td>29.50</td><td>31.98</td></tr><tr><td>(0.02,0.02,720)</td><td>71.16</td><td>71.60</td><td>71.60</td><td>71.74</td><td>71.83</td><td>71.83</td><td>72.79</td></tr><tr><td>(0.04,0.04,720)</td><td>46.32</td><td>46.41</td><td>46.50</td><td>46.53</td><td>46.52</td><td>46.52</td><td>48.25</td></tr><tr><td>(0.06,0.06,720)</td><td>17.50</td><td>17.11</td><td>17.10</td><td>17.04</td><td>17.01</td><td>16.99</td><td>17.75</td></tr><tr><td>(0.02,0.02,1800)</td><td>27.98</td><td>28.38</td><td>28.41</td><td>28.52</td><td>28.59</td><td>28.61</td><td>29.06</td></tr><tr><td>(0.04,0.04,1800)</td><td>17.02</td><td>17.17</td><td>17.22</td><td>17.25</td><td>17.27</td><td>17.28</td><td>17.12</td></tr><tr><td>(0.06,0.06,1800)</td><td>6.160</td><td>6.036</td><td>6.036</td><td>6.016</td><td>6.012</td><td>6.005</td><td>5.76</td></tr></tbody></table>

The COMSOL numerical solution is also analyzed in terms of the dimensionless temperature, $\theta=(T_{in}-T)/T_{in}$, relative deviations with respect to the converged GITT results, Table V. The deviations are quite satisfactory, being throughout most of the time and space domains below 5%.
<table><caption>Table V – Dimensionless temperature (COMSOL) relative deviations with respect to GITT solution with N=M=600 at different positions.</caption><thead><tr><th>M<br>(x, y, t)</th><th>GITT<br>(M=N=600)</th><th>COMSOL<br>(Extra Fine)</th><th>Relative deviation<br>(%)</th></tr></thead><tbody><tr><td>(0.02,0.02,360)</td><td>0.07359</td><td>0.07969</td><td>8.29</td></tr><tr><td>(0.04,0.04,360)</td><td>0.2934</td><td>0.2766</td><td>5.71</td></tr><tr><td>(0.06,0.06,360)</td><td>0.7050</td><td>0.6802</td><td>3.51</td></tr><tr><td>(0.02,0.02,720)</td><td>0.2817</td><td>0.2721</td><td>3.40</td></tr><tr><td>(0.04,0.04,720)</td><td>0.5348</td><td>0.5175</td><td>3.24</td></tr><tr><td>(0.06,0.06,720)</td><td>0.8301</td><td>0.8225</td><td>0.920</td></tr><tr><td>(0.02,0.02,1080)</td><td>0.4668</td><td>0.4534</td><td>2.86</td></tr><tr><td>(0.04,0.04,1080)</td><td>0.6716</td><td>0.6659</td><td>0.839</td></tr><tr><td>(0.06,0.06,1080)</td><td>0.8842</td><td>0.8850</td><td>0.0875</td></tr><tr><td>(0.02,0.02,1800)</td><td>0.7139</td><td>0.7094</td><td>0.638</td></tr><tr><td>(0.04,0.04,1800)</td><td>0.8272</td><td>0.8288</td><td>0.187</td></tr><tr><td>(0.06,0.06,1800)</td><td>0.9399</td><td>0.9424</td><td>0.261</td></tr></tbody></table>

Figures 2.a-c further illustrate the COMSOL [34] mesh refinement analysis, according to the selected options shown in Table III, plotting the time evolution of the dimensionless temperature, $\theta=(T_{in}-T)/T_{in}$, relative deviations with respect to the present GITT results converged to the fourth digit, with M=N=600 terms, at three different positions. Again, one may observe that the agreement is quite satisfactory for the different meshes, with relative deviations falling below 5% along most of the time evolution. It is also observable, for each spatial position considered, the time interval when gradients and time rates are locally more significant, and the numerical approximation leads to a relative deviation increase.
<table><caption>Table V – Dimensionless temperature (COMSOL) relative deviations with respect to GITT solution with $N=M=600$ at different positions.</caption><thead><tr><th>M<br>(x, y, t)</th><th>GITT<br>(M=N=600)</th><th>COMSOL<br>(Extra Fine)</th><th>Relative deviation<br>(%)</th></tr></thead><tbody><tr><td>(0.02,0.02,360)</td><td>0.07359</td><td>0.07969</td><td>8.29</td></tr><tr><td>(0.04,0.04,360)</td><td>0.2934</td><td>0.2766</td><td>5.71</td></tr><tr><td>(0.06,0.06,360)</td><td>0.7050</td><td>0.6802</td><td>3.51</td></tr><tr><td>(0.02,0.02,720)</td><td>0.2817</td><td>0.2721</td><td>3.40</td></tr><tr><td>(0.04,0.04,720)</td><td>0.5348</td><td>0.5175</td><td>3.24</td></tr><tr><td>(0.06,0.06,720)</td><td>0.8301</td><td>0.8225</td><td>0.920</td></tr><tr><td>(0.02,0.02,1080)</td><td>0.4668</td><td>0.4534</td><td>2.86</td></tr><tr><td>(0.04,0.04,1080)</td><td>0.6716</td><td>0.6659</td><td>0.839</td></tr><tr><td>(0.06,0.06,1080)</td><td>0.8842</td><td>0.8850</td><td>0.0875</td></tr><tr><td>(0.02,0.02,1800)</td><td>0.7139</td><td>0.7094</td><td>0.638</td></tr><tr><td>(0.04,0.04,1800)</td><td>0.8272</td><td>0.8288</td><td>0.187</td></tr><tr><td>(0.06,0.06,1800)</td><td>0.9399</td><td>0.9424</td><td>0.261</td></tr></tbody></table>

Figures 2.a-c further illustrate the COMSOL [34] mesh refinement analysis, according to the selected options shown in Table III, plotting the time evolution of the dimensionless temperature, $\theta=(T_{in}-T)/T_{in}$, relative deviations with respect to the present GITT results converged to the fourth digit, with $M=N=600$ terms, at three different positions. Again, one may observe that the agreement is quite satisfactory for the different meshes, with relative deviations falling below 5% along most of the time evolution. It is also observable, for each spatial position considered, the time interval when gradients and time rates are locally more significant, and the numerical approximation leads to a relative deviation increase.
<figure><img src="image_3.png" /></figure>

Figures 2a-c – COMSOL [34] dimensionless temperature relative deviations for different mesh refinement options, Table III, with respect to converged GITT results at positions (x,y): (a) (0.02,0.02); (b) (0.04,0.04); (c) (0.06,0.06).
<figure><img src="image_3.png" /></figure>

Figures 2a-c – COMSOL [34] dimensionless temperature relative deviations for different mesh refinement options, Table III, with respect to converged GITT results at positions (x,y): (a) (0.02,0.02); (b) (0.04,0.04); (c) (0.06,0.06).
Figure 3 provides temperature profiles at different time values along the y variable for the vertical line with x=0.04, which crosses both the substrate and the anisotropic inclusion right at the midline, thus capturing the transitions between the two materials. The present converged GITT results are shown in solid lines, in different colours for each time value (t=180, 360, 720, and 1800 s), while the COMSOL results with the Extra Fine mesh refinement option are provided in dotted lines. Again, the excellent agreement between the numerical and hybrid numerical-analytical solutions is evident, with practically full adherence to the graph scale. The central inclusion thermal behaviour, due to its larger anisotropic conductivity, is noticeable in terms of the reduced temperature gradients in this region.

<figure><img src="image_4.png" /><figcaption>Figure 3 – Temperature profiles along y-variable for vertical line at x=a/2=0.04m (solid lines – GITT with M=600; dashed lines – COMSOL with Extra Fine mesh)</figcaption></figure>

Next, we consider two asymmetric situations of a rectangular central inclusion, with anisotropic thermal conductivity. The first asymmetric situation is a rectangular insert in the middle of the substrate plate, with sides $l_x=0.015m$ and $l_y=0.025m$, and thermal conductivity components $k_{11}=k_{22}=26 W/m^°C$, $k_{12}=k_{21}=8.667 W/m^°C$, as in the previous square insert case. The second asymmetric situation is the same rectangular insert in size but with different thermal conductivity components, $k_{11}=13 W/m^°C$, $k_{22}=52 W/m^°C$, $k_{12}=k_{21}=8.667 W/m^°C$. The remaining thermophysical properties assume the same values as for the previous square central region. Figures 4a,b provide the temperature profiles at different time values along the x and y variables for the horizontal line with y=0.04 m the vertical line with x=0.04 m, respectively, which crosses both the substrate and the anisotropic inclusion right at the midline, thus capturing the transitions
Figure 3 provides temperature profiles at different time values along the y variable for the vertical line with x=0.04, which crosses both the substrate and the anisotropic inclusion right at the midline, thus capturing the transitions between the two materials. The present converged GITT results are shown in solid lines, in different colours for each time value (t=180, 360, 720, and 1800 s), while the COMSOL results with the Extra Fine mesh refinement option are provided in dotted lines. Again, the excellent agreement between the numerical and hybrid numerical-analytical solutions is evident, with practically full adherence to the graph scale. The central inclusion thermal behaviour, due to its larger anisotropic conductivity, is noticeable in terms of the reduced temperature gradients in this region.

<figure><img src="image_4.png" /><figcaption>Figure 3 – Temperature profiles along y-variable for vertical line at x=a/2=0.04m (solid lines – GITT with M=600; dashed lines – COMSOL with Extra Fine mesh)</figcaption></figure>

Next, we consider two asymmetric situations of a rectangular central inclusion, with anisotropic thermal conductivity. The first asymmetric situation is a rectangular insert in the middle of the substrate plate, with sides $l_x=0.015m$ and $l_y=0.025m$, and thermal conductivity components $k_{11}=k_{22}=26 W/m^°C$, $k_{12}=k_{21}=8.667 W/m^°C$, as in the previous square insert case. The second asymmetric situation is the same rectangular insert in size but with different thermal conductivity components, $k_{11}=13 W/m^°C$, $k_{22}=52 W/m^°C$, $k_{12}=k_{21}=8.667 W/m^°C$. The remaining thermophysical properties assume the same values as for the previous square central region. Figures 4a,b provide the temperature profiles at different time values along the x and y variables for the horizontal line with y=0.04 m the vertical line with x=0.04 m, respectively, which crosses both the substrate and the anisotropic inclusion right at the midline, thus capturing the transitions
between the two materials. Figures 5a,b provide again the temperature profiles along the same horizontal and vertical lines, but now for the other case with different diagonal thermal conductivity values. The GITT results with M=500 terms are shown in solid lines, in different colours for each time value (t=180, 360, 720, and 1800 s), while the COMSOL results with the Extra Fine mesh refinement option are provided in dotted lines, again with excellent agreement to the graph scale with the integral transform benchmark results.

<figure><img src="image_5.png" /><figcaption>Figure 4 – Temperature profiles along (a) horizontal line at y=b/2=0.04m and (b) vertical line at x=a/2=0.04m with k₁₁=k₂₂=26 W/m·°C, k₁₂= k₂₁=8.667 W/m·°C (solid lines: GITT M=500; dotted lines: COMSOL Extra Fine mesh)</figcaption></figure>
between the two materials. Figures 5a,b provide again the temperature profiles along the same horizontal and vertical lines, but now for the other case with different diagonal thermal conductivity values. The GITT results with M=500 terms are shown in solid lines, in different colours for each time value (t=180, 360, 720, and 1800 s), while the COMSOL results with the Extra Fine mesh refinement option are provided in dotted lines, again with excellent agreement to the graph scale with the integral transform benchmark results.

<figure><img src="image_5.png" /><figcaption>Figure 4 – Temperature profiles along (a) horizontal line at y=b/2=0.04m and (b) vertical line at x=a/2=0.04m with k₁₁=k₂₂=26 W/m·°C, k₁₂= k₂₁=8.667 W/m·°C (solid lines: GITT M=500; dotted lines: COMSOL Extra Fine mesh)</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5 – Temperature profiles along (a) horizontal line at y=b/2=0.04m and (b) vertical line at x=a/2=0.04m with k₁₁=13 W/m°C, k₂₂=52 W/m°C, k₁₂= k₂₁=8.667 W/m°C (solid lines – GITT M=500; dotted lines – COMSOL Extra Fine mesh)</figcaption></figure>

Finally, Figures 6a,b provide the isotherms within the two regions, substrate and central inclusion, at *t*=900 seconds, as obtained from the GITT solution with M=600 terms, for the squared and rectangular inserts with different thermal conductivities. From the dashed regions one can observe the fairly uniform temperature distribution within the central inclusions, with the corresponding deformations due to the anisotropic behavior.
<figure><img src="image_6.png" /><figcaption>Figure 5 – Temperature profiles along (a) horizontal line at y=b/2=0.04m and (b) vertical line at x=a/2=0.04m with k₁₁=13 W/m°C, k₂₂=52 W/m°C, k₁₂= k₂₁=8.667 W/m°C (solid lines – GITT M=500; dotted lines – COMSOL Extra Fine mesh)</figcaption></figure>

Finally, Figures 6a,b provide the isotherms within the two regions, substrate and central inclusion, at *t*=900 seconds, as obtained from the GITT solution with M=600 terms, for the squared and rectangular inserts with different thermal conductivities. From the dashed regions one can observe the fairly uniform temperature distribution within the central inclusions, with the corresponding deformations due to the anisotropic behavior.
<figure><img src="image_7.png" /><figcaption>Figure 6 – Isotherms of GITT solution (M=600) at t=900 s. Dashed line is the central anisotropic inclusion: (a)square, k₁₁=k₂₂=26 W/m°C, k₁₂=k₂₁=8.667 W/m°C; (b)rectangle, k₁₁=13 W/m°C, k₂₂=52 W/m°C, k₁₂=k₂₁=8.667 W/m°C</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6 – Isotherms of GITT solution (M=600) at t=900 s. Dashed line is the central anisotropic inclusion: (a)square, k₁₁=k₂₂=26 W/m°C, k₁₂=k₂₁=8.667 W/m°C; (b)rectangle, k₁₁=13 W/m°C, k₂₂=52 W/m°C, k₁₂=k₂₁=8.667 W/m°C</figcaption></figure>
CONCLUSIONS

The Generalized Integral Transform Technique (GITT) was further advanced to obtain hybrid numerical-analytical solutions for transient heat conduction in multidimensional heterogeneous anisotropic domains. Based on previous works, a biorthogonal eigenfunction expansion is proposed that carries along the information on the anisotropic diffusive terms. As for the heterogeneous behavior, a single domain reformulation strategy rewrites the problem with multiple regions or materials as one single energy equation with space variable coefficients and source terms. Besides, the GITT itself is proposed for the solution of the biorthogonal eigenvalue problems pair with the heterogeneity incorporated into the coefficients, thus transforming the differential eigenvalue problems into matrix eigensystem analysis. Simpler biorthogonal or self-adjoint auxiliary eigenvalue problems are considered in writing the expansions for the original unknown eigenfunctions. A test case consisting of an isotropic rectangular substrate with a central anisotropic inclusion is considered more closely and used to illustrate the excellent convergence behavior of the temperature field based on the biorthogonal eigenvalue problems. Besides, this error-controlled solution is employed as benchmark reference results to inspect the mesh refinement improvement in finite element solutions for the same test-case, as obtained from the well-known COMSOL Multiphysics code [34]. The proposed approach is directly applicable to handling heat conduction in functionally graded heterogeneous media and thermal metamaterials modelling. Besides, it opens avenues of extension possibilities in the analysis of heat transfer in irregular anisotropic domains [36-37] and in the identification of thermophysical properties in heterogeneous anisotropic materials [38-39].

ACKNOWLEDGEMENTS

The authors wish to express their gratitude for the partial financial support from
the sponsoring agencies CAPES (PROCAD-Defesa), CNPq, and FAPERJ, all of them in
Brazil. ML would like to acknowledge the hospitality of the Mechanical Eng. Dept. of
COPPE/UFRJ and the financial support of FAPERJ and the French Embassy in Brazil,
during his two months visit to UFRJ in 2022. AEHB would like to acknowledge the
CONCLUSIONS

The Generalized Integral Transform Technique (GITT) was further advanced to obtain hybrid numerical-analytical solutions for transient heat conduction in multidimensional heterogeneous anisotropic domains. Based on previous works, a biorthogonal eigenfunction expansion is proposed that carries along the information on the anisotropic diffusive terms. As for the heterogeneous behavior, a single domain reformulation strategy rewrites the problem with multiple regions or materials as one single energy equation with space variable coefficients and source terms. Besides, the GITT itself is proposed for the solution of the biorthogonal eigenvalue problems pair with the heterogeneity incorporated into the coefficients, thus transforming the differential eigenvalue problems into matrix eigensystem analysis. Simpler biorthogonal or self-adjoint auxiliary eigenvalue problems are considered in writing the expansions for the original unknown eigenfunctions. A test case consisting of an isotropic rectangular substrate with a central anisotropic inclusion is considered more closely and used to illustrate the excellent convergence behavior of the temperature field based on the biorthogonal eigenvalue problems. Besides, this error-controlled solution is employed as benchmark reference results to inspect the mesh refinement improvement in finite element solutions for the same test-case, as obtained from the well-known COMSOL Multiphysics code [34]. The proposed approach is directly applicable to handling heat conduction in functionally graded heterogeneous media and thermal metamaterials modelling. Besides, it opens avenues of extension possibilities in the analysis of heat transfer in irregular anisotropic domains [36-37] and in the identification of thermophysical properties in heterogeneous anisotropic materials [38-39].

ACKNOWLEDGEMENTS

The authors wish to express their gratitude for the partial financial support from
the sponsoring agencies CAPES (PROCAD-Defesa), CNPq, and FAPERJ, all of them in
Brazil. ML would like to acknowledge the hospitality of the Mechanical Eng. Dept. of
COPPE/UFRJ and the financial support of FAPERJ and the French Embassy in Brazil,
during his two months visit to UFRJ in 2022. AEHB would like to acknowledge the
hospitality and the financial support of URCA during his three months stay at the
Université de Reims in 2023.

DEDICATION

This work is dedicated to the memory of Prof. Darrell W. Pepper, co-founder and
editor-in-chief of the Computational Thermal Sciences journal. Through his outstanding
perception of the computational engineering sciences field, his continuous efforts in
benchmarking numerical computations, and his appreciation of the importance and merits
of the meshless methods alternative, he has always been a sincere supporter of our work
on the integral transform method, itself a benchmarking tool and a meshless method,
throughout the years.

REFERENCES

1. M. N. Ozisik, Heat Conduction, John Wiley and Sons, New York (1980).

2. G.P. Mulholland, B.P. Gupta, Heat transfer in a three-dimensional anisotropic solid of arbitrary shape, J. Heat Transfer 99 (1977) 135–137.

3. Y.P. Chang, Analytical solution for heat conduction in anisotropic media in infinite semi-infinite, and two-place bounded regions, Int. J. Heat Mass Transfer 20 (1977) 1019–1028.

4. K.C. Poon, R.C.H. Tsou, Y.P. Chang, Solution of anisotropic problems of first class by coordinate-transformation, J. Heat Transfer 101 (1979) 340–345.

5. Sladek, J.; Sladek, V.; Hellmich, Ch.; Eberhardsteiner, J. (2007a): Heat conduction analysis of 3D axisymmetric and anisotropic FGM bodies by meshless local Petrov-Galerkin method, Computational Mechanics, 39: 323-333.

6. M.I. Azis and D.L. Clements, A boundary element method for transient heat conduction problem of nonhomogeneous anisotropic materials. Far East J Math Sci 2014;89(1):51–67.

7. Robert Schittny, Muamer Kadic, Sebastien Guenneau, and Martin Wegener. Experiments on transformation thermodynamics: Molding the flow of heat. Physical Review Letters, 110(19), 2013.
hospitality and the financial support of URCA during his three months stay at the
Université de Reims in 2023.

DEDICATION

This work is dedicated to the memory of Prof. Darrell W. Pepper, co-founder and
editor-in-chief of the Computational Thermal Sciences journal. Through his outstanding
perception of the computational engineering sciences field, his continuous efforts in
benchmarking numerical computations, and his appreciation of the importance and merits
of the meshless methods alternative, he has always been a sincere supporter of our work
on the integral transform method, itself a benchmarking tool and a meshless method,
throughout the years.

REFERENCES

1. M. N. Ozisik, Heat Conduction, John Wiley and Sons, New York (1980).

2. G.P. Mulholland, B.P. Gupta, Heat transfer in a three-dimensional anisotropic solid of arbitrary shape, J. Heat Transfer 99 (1977) 135–137.

3. Y.P. Chang, Analytical solution for heat conduction in anisotropic media in infinite semi-infinite, and two-place bounded regions, Int. J. Heat Mass Transfer 20 (1977) 1019–1028.

4. K.C. Poon, R.C.H. Tsou, Y.P. Chang, Solution of anisotropic problems of first class by coordinate-transformation, J. Heat Transfer 101 (1979) 340–345.

5. Sladek, J.; Sladek, V.; Hellmich, Ch.; Eberhardsteiner, J. (2007a): Heat conduction analysis of 3D axisymmetric and anisotropic FGM bodies by meshless local Petrov-Galerkin method, Computational Mechanics, 39: 323-333.

6. M.I. Azis and D.L. Clements, A boundary element method for transient heat conduction problem of nonhomogeneous anisotropic materials. Far East J Math Sci 2014;89(1):51–67.

7. Robert Schittny, Muamer Kadic, Sebastien Guenneau, and Martin Wegener. Experiments on transformation thermodynamics: Molding the flow of heat. Physical Review Letters, 110(19), 2013.
8. Supradeep Narayana, Salvatore Savo, and Yuki Sato, Transient heat flux shielding using thermal metamaterials, Appl. Phys. Lett. 102, 201904 (2013).

9. J.P. Huang, Theoretical Thermotics, Springer, Singapore, 2020.

10. Muhammad Imra, Liangchi Zhang, Asit Kumar Gain, Advanced thermal metamaterial design for temperature control at the cloaked region, Nature Scientific Reports, (2020) 10:11763

11. M.D. Mikhailov and M.N. Ozisik, Unified Analysis and Solutions of Heat and Mass Diffusion, John Wiley, New York (1984).

12. M.D. Mikhailov and M.N. Ozisik, On general solution of heat conduction in an anisotropic medium. Letters in Heat and Mass Transfer, 8, 329-335 (1981).

13. M.D. Mikhailov and M.N. Ozisik, Unified finite integral transform analysis, J. Franklin Institute, 321, 299-307 (1986).

14. Cotta, R.M. (1990) Hybrid Numerical-Analytical Approach to Nonlinear Diffusion Problems. Num. Heat Transfer, Part B, 127: 217-226.

15. Serfaty, R., Cotta, R.M., Integral transform solutions of diffusion problems with nonlinear equation coefficients, Int. Comm. Heat & Mass Transfer, 17 (1990) 851-864.

16. Cotta, R.M. (1993) Integral Transforms in Computational Heat and Fluid Flow. CRC Press, Boca Raton, FL.

17. Cotta, R.M. (1994) Benchmark Results in Computational Heat and Fluid Flow: - The Integral Transform Method. Int. J. Heat Mass Transfer (Invited Paper), 37:381-394.

18. Cotta, R.M., Mikhailov, M.D. (1997) Heat Conduction: Lumped Analysis, Integral Transforms, Symbolic Computation. Wiley, Chichester.

19. Mikhailov, M.D. and R.M. Cotta, Unified Integral Transform Method, Journal of the Brazilian Society of Mechanical Sciences and Eng. (Invited Paper), Vol. 12, No.3, pp. 301-310, 1990.

20. Cotta, R.M., Mikhailov, M.D. (2006) Hybrid Methods and Symbolic Computations. In: Handbook of Numerical Heat Transfer, 2nd edition, Chapter 16, Eds. W.J. Minkowycz, E.M. Sparrow, and J.Y. Murthy, John Wiley, New York.

21. Cotta, R.M., D.C. Knupp, C. P. Naveira-Cotta, L.A. Sphaier, and J.N.N. Quaresma, "Unified Integral Transforms Algorithm for Solving Multidimensional Nonlinear Convection-Diffusion Problems", Num. Heat Transfer, part A - Applications, V.63, pp.1-27, 2013.
8. Supradeep Narayana, Salvatore Savo, and Yuki Sato, Transient heat flux shielding using thermal metamaterials, Appl. Phys. Lett. 102, 201904 (2013).

9. J.P. Huang, Theoretical Thermotics, Springer, Singapore, 2020.

10. Muhammad Imra, Liangchi Zhang, Asit Kumar Gain, Advanced thermal metamaterial design for temperature control at the cloaked region, Nature Scientific Reports, (2020) 10:11763

11. M.D. Mikhailov and M.N. Ozisik, Unified Analysis and Solutions of Heat and Mass Diffusion, John Wiley, New York (1984).

12. M.D. Mikhailov and M.N. Ozisik, On general solution of heat conduction in an anisotropic medium. Letters in Heat and Mass Transfer, 8, 329-335 (1981).

13. M.D. Mikhailov and M.N. Ozisik, Unified finite integral transform analysis, J. Franklin Institute, 321, 299-307 (1986).

14. Cotta, R.M. (1990) Hybrid Numerical-Analytical Approach to Nonlinear Diffusion Problems. Num. Heat Transfer, Part B, 127: 217-226.

15. Serfaty, R., Cotta, R.M., Integral transform solutions of diffusion problems with nonlinear equation coefficients, Int. Comm. Heat & Mass Transfer, 17 (1990) 851-864.

16. Cotta, R.M. (1993) Integral Transforms in Computational Heat and Fluid Flow. CRC Press, Boca Raton, FL.

17. Cotta, R.M. (1994) Benchmark Results in Computational Heat and Fluid Flow: - The Integral Transform Method. Int. J. Heat Mass Transfer (Invited Paper), 37:381-394.

18. Cotta, R.M., Mikhailov, M.D. (1997) Heat Conduction: Lumped Analysis, Integral Transforms, Symbolic Computation. Wiley, Chichester.

19. Mikhailov, M.D. and R.M. Cotta, Unified Integral Transform Method, Journal of the Brazilian Society of Mechanical Sciences and Eng. (Invited Paper), Vol. 12, No.3, pp. 301-310, 1990.

20. Cotta, R.M., Mikhailov, M.D. (2006) Hybrid Methods and Symbolic Computations. In: Handbook of Numerical Heat Transfer, 2nd edition, Chapter 16, Eds. W.J. Minkowycz, E.M. Sparrow, and J.Y. Murthy, John Wiley, New York.

21. Cotta, R.M., D.C. Knupp, C. P. Naveira-Cotta, L.A. Sphaier, and J.N.N. Quaresma, "Unified Integral Transforms Algorithm for Solving Multidimensional Nonlinear Convection-Diffusion Problems", Num. Heat Transfer, part A - Applications, V.63, pp.1-27, 2013.
22. R.M. Cotta, D.C. Knupp, and C.P. Naveira-Cotta, Analytical Heat and Fluid Flow in Microchannels and Microsystems, Mechanical Eng. Series, Springer-Verlag, 2016.

23. Cotta, R.M., Knupp, D.C., Quaresma, J.N.N. (2018) Analytical Methods in Heat Transfer. In: Handbook of Thermal Science and Engineering, F.A. Kulacki et al., Eds., Springer.

24. Cotta, R.M., Naveira-Cotta, C.P., Knupp, D.C., Zotin, J.L.Z., Pontes, P.C., Almeida, A.P. (2018) Recent Advances in Computational-Analytical Integral Transforms for Convection-Diffusion Problems. Heat & Mass Transfer, Invited Paper, 54:2475–2496.

25. Mikhailov, M.D., Cotta, R.M., Integral transform method for eigenvalue problems, Comm. Num. Meth. Eng., 10 (1994) 827-835.

26. Naveira-Cotta, C.P., R.M. Cotta, H.R.B. Orlande, and O. Fudym, "Eigenfunction Expansions for Transient Diffusion in Heterogeneous Media", Int. J. Heat and Mass Transfer, V.52, pp.5029-5039, 2009.

27. Knupp, D.C., C.P. Naveira-Cotta, J.V.C. Ayres, R.M. Cotta, and H.R.B. Orlande, "Theoretical-Experimental Analysis of Heat Transfer in Nonhomogeneous Solids Via Improved Lumped Formulation, Integral Transforms and Infrared Thermography", Int. J. Thermal Sciences, V.62, pp.71-84, 2012.

28. Knupp, D.C., Naveira-Cotta, C.P., Cotta, R.M. (2012) Theoretical Analysis of Conjugated Heat Transfer with a Single Domain Formulation and Integral Transforms. Int. Comm. Heat & Mass Transfer, V.39, no.3, pp.355-362.

29. Knupp, D.C., Naveira-Cotta, C.P., Cotta, R.M. (2013) Conjugated convection-conduction analysis in micro-channels with axial diffusion effects and a single domain formulation, ASME Journal of Heat Transfer, Special Issue on Micro/Nanoscale Heat and Mass Transfer. V.135, 091401-091411.

30. Knupp, D.C., Cotta, R.M., Naveira-Cotta, C.P., Kakaç, S. (2015) Transient conjugated heat transfer in microchannels: Integral transforms with single domain formulation. Int. J. Thermal Sciences, Vol. 88, pp. 248-257.

31. Knupp, D.C., Cotta, R.M., Naveira-Cotta, C.P. (2015) Fluid flow and conjugated heat transfer in arbitrarily shaped channels via single domain formulation and integral transforms. Int. J. Heat Mass Transfer, Vol. 82, pp. 479-489.

32. Lisboa, K.M., Su, J., Cotta, R.M. (2018) Single Domain Integral Transforms Analysis of Natural Convection in Cavities Partially Filled with Heat Generating Porous Medium, Num. Heat Transfer, Part A – Applications, V.74, no.3, pp.1068-1086.
22. R.M. Cotta, D.C. Knupp, and C.P. Naveira-Cotta, Analytical Heat and Fluid Flow in Microchannels and Microsystems, Mechanical Eng. Series, Springer-Verlag, 2016.

23. Cotta, R.M., Knupp, D.C., Quaresma, J.N.N. (2018) Analytical Methods in Heat Transfer. In: Handbook of Thermal Science and Engineering, F.A. Kulacki et al., Eds., Springer.

24. Cotta, R.M., Naveira-Cotta, C.P., Knupp, D.C., Zotin, J.L.Z., Pontes, P.C., Almeida, A.P. (2018) Recent Advances in Computational-Analytical Integral Transforms for Convection-Diffusion Problems. Heat & Mass Transfer, Invited Paper, 54:2475–2496.

25. Mikhailov, M.D., Cotta, R.M., Integral transform method for eigenvalue problems, Comm. Num. Meth. Eng., 10 (1994) 827-835.

26. Naveira-Cotta, C.P., R.M. Cotta, H.R.B. Orlande, and O. Fudym, "Eigenfunction Expansions for Transient Diffusion in Heterogeneous Media", Int. J. Heat and Mass Transfer, V.52, pp.5029-5039, 2009.

27. Knupp, D.C., C.P. Naveira-Cotta, J.V.C. Ayres, R.M. Cotta, and H.R.B. Orlande, "Theoretical-Experimental Analysis of Heat Transfer in Nonhomogeneous Solids Via Improved Lumped Formulation, Integral Transforms and Infrared Thermography", Int. J. Thermal Sciences, V.62, pp.71-84, 2012.

28. Knupp, D.C., Naveira-Cotta, C.P., Cotta, R.M. (2012) Theoretical Analysis of Conjugated Heat Transfer with a Single Domain Formulation and Integral Transforms. Int. Comm. Heat & Mass Transfer, V.39, no.3, pp.355-362.

29. Knupp, D.C., Naveira-Cotta, C.P., Cotta, R.M. (2013) Conjugated convection-conduction analysis in micro-channels with axial diffusion effects and a single domain formulation, ASME Journal of Heat Transfer, Special Issue on Micro/Nanoscale Heat and Mass Transfer. V.135, 091401-091411.

30. Knupp, D.C., Cotta, R.M., Naveira-Cotta, C.P., Kakaç, S. (2015) Transient conjugated heat transfer in microchannels: Integral transforms with single domain formulation. Int. J. Thermal Sciences, Vol. 88, pp. 248-257.

31. Knupp, D.C., Cotta, R.M., Naveira-Cotta, C.P. (2015) Fluid flow and conjugated heat transfer in arbitrarily shaped channels via single domain formulation and integral transforms. Int. J. Heat Mass Transfer, Vol. 82, pp. 479-489.

32. Lisboa, K.M., Su, J., Cotta, R.M. (2018) Single Domain Integral Transforms Analysis of Natural Convection in Cavities Partially Filled with Heat Generating Porous Medium, Num. Heat Transfer, Part A – Applications, V.74, no.3, pp.1068-1086.
33. Almeida, A.P., C.P. Naveira-Cotta, and R.M. Cotta, "Transient Multidimensional Heat Conduction in Heterogeneous Media: Integral Transform with Single Domain Formulation", Int. Comm. Heat & Mass Transfer, V.117, pp.104792 (13p.), 2020.

34. COMSOL Multiphysics Reference Manual, version 5.6, COMSOL Inc, 2020.
www.comsol.com

35. Wolfram, S., Mathematica, version 13.2, Wolfram Research Inc., Champaign, IL, 2022.

36. Aparecido, J.B., Cotta, R.M., Ozisik, M.N. (1989) Analytical Solutions to Two-Dimensional Diffusion Type Problems in Irregular Geometries, J. of the Franklin Institute, V. 326, pp. 421-434.

37. Sphaier, L. A., Cotta, R.M. (2022) Analytical and Hybrid Solutions of Diffusion Problems within Arbitrarily Shaped Regions via Integral Transforms", Computational Mechanics, V.29, no.3, pp.265-276.

38. Naveira-Cotta, C.P., Cotta, R.M., and Orlande, H.R.B. (2011) Inverse analysis with integral transformed temperature fields for identification of thermophysical properties functions in heterogeneous media, International Journal of Heat and Mass Transfer. V.54, pp.1506-1519.

39. Knupp, D.C., Naveira-Cotta, C.P., Ayres, J.V.C, Cotta, R.M., and Orlande, H.R.B. (2012) Theoretical-experimental analysis of heat transfer in non-homogeneous solids via improved lumped formulation, integral transforms, and infrared thermography, International Journal of Thermal Sciences. V.62, pp.71-84.
33. Almeida, A.P., C.P. Naveira-Cotta, and R.M. Cotta, "Transient Multidimensional Heat Conduction in Heterogeneous Media: Integral Transform with Single Domain Formulation", Int. Comm. Heat & Mass Transfer, V.117, pp.104792 (13p.), 2020.

34. COMSOL Multiphysics Reference Manual, version 5.6, COMSOL Inc, 2020.
www.comsol.com

35. Wolfram, S., Mathematica, version 13.2, Wolfram Research Inc., Champaign, IL, 2022.

36. Aparecido, J.B., Cotta, R.M., Ozisik, M.N. (1989) Analytical Solutions to Two-Dimensional Diffusion Type Problems in Irregular Geometries, J. of the Franklin Institute, V. 326, pp. 421-434.

37. Sphaier, L. A., Cotta, R.M. (2022) Analytical and Hybrid Solutions of Diffusion Problems within Arbitrarily Shaped Regions via Integral Transforms", Computational Mechanics, V.29, no.3, pp.265-276.

38. Naveira-Cotta, C.P., Cotta, R.M., and Orlande, H.R.B. (2011) Inverse analysis with integral transformed temperature fields for identification of thermophysical properties functions in heterogeneous media, International Journal of Heat and Mass Transfer. V.54, pp.1506-1519.

39. Knupp, D.C., Naveira-Cotta, C.P., Ayres, J.V.C, Cotta, R.M., and Orlande, H.R.B. (2012) Theoretical-experimental analysis of heat transfer in non-homogeneous solids via improved lumped formulation, integral transforms, and infrared thermography, International Journal of Thermal Sciences. V.62, pp.71-84.