# Simulation of effective scale-size dependent heat conduction in rigid microgeometries

Mario Setta, Eddie Wadbro, Grigor Nika

► To cite this version:

Mario Setta, Eddie Wadbro, Grigor Nika. Simulation of effective scale-size dependent heat conduction in rigid microgeometries. 2025. hal-05299481

HAL Id: hal-05299481

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

Preprint submitted on 6 Oct 2025

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
Simulation of effective scale-size
dependent heat conduction in rigid microgeometries

Mario Settaª, Eddie Wadbroª, Grigor Nikaª,*

ªDepartment of Mathematics and Computer Science, Karlstad University, SE-651 88, Karlstad, Sweden

**Abstract**

We present homogenization and simulation results for an enhanced heat equation model that captures thermal scale-size effects through higher-gradient corrections involving characteristic internal lengths. The resulting equation is a fourth-order parabolic equation that incorporates thermal scale effects inherent to microstructured materials. We derive effective thermal coefficients for the time-stationary problem using asymptotic homogenization. This enables accurate simulation via a quadratic B-spline-based finite element approach. Our results quantify the influence of microstructure shape and volume fraction on the effective thermal behavior, demonstrating how scale-size-induced phenomena critically affect heat transport in micro- and nanoscale devices.

Keywords: scale-size thermal effects, homogenization, quadratic splines, gradient of entropy, thermodynamics

MSC 2020: 80A17, 80M35, 80M40, 74F05, 74M25

1. Introduction

Over the last 30 years, the trend towards smaller electronic components and integrated circuits on chip
devices has attracted significant attention and been a focus of research and development. Miniaturization
has led to compact, high-performance devices such as smartphones and OLED screens, the recent *boom* in
artificial intelligence, robotics, and many other fields. Simultaneously, as improved additive manufacturing
techniques enabled the aforementioned reduction in the length scale of components of devices, new challenges
emerged as the size of the components or the dimension of the device approached the molecular level. At
such small scales (see Fig. 1 for a perspective on the different scales involved), classical continuum mechanics
loses its accuracy as its core assumption of scale-invariant material properties breaks down and significant
scale size effects emerge that alter the behavior of the device [2, 5, 13, 21, 22, 40, 25, 24]. Due to the high
computational cost of atomistic simulations, generalized continuum theories [29, 30, 44, 45] incorporating
higher-gradient corrections depending of characteristic internal lengths offer a bridge connecting continuum
and atomistic physics [23]. The application of such gradient-based corrections has an established history in
the field of rarefied gas dynamics in modeling convective heat transfer in microdomains (with characteristic
lengths ranging from 0.1 to 100 µm). Recent studies have further explored coupling fluid dynamics with
higher-gradient heat equations (see, e.g., [17, 43]). In this *slip flow* regime, where classical Fourier conduction
is inadequate, the key parameter is the Knudsen number, Kn. Defined as the ratio of the mean free path to
a characteristic length scale (see, e.g., [4], [8]), Kn governs the transport behavior as non-continuum effects
become significant with increasing Kn. This dimensionless number is directly related to the internal length
scales that describe these thermal scale-size effects.

<footnote>*Corresponding author</footnote>

<footnote>Email addresses: mario.setta@kau.se (Mario Setta), eddie.wadbro@kau.se (Eddie Wadbro), grigor.nika@kau.se
(Grigor Nika)</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: The characteristic length scale of the microstructure of the bespoken devices in perspective. The image is a modified version of the image in [8].</figcaption></figure>

This work aims to (i) derive a homogenized model for the *gradient of entropy* heat equation taking into account the size of the representative volume element; (ii) provide a numerical modeling framework for the simulation of the effective coefficients based on quadratic B-splines in order to account for the $C^1$ interpolation between the elements that is required due to the higher-gradient terms; and (iii) perform numerical simulations for different size representative volume elements with different volume fractions and quantify their effect on the homogenized coefficients.

We derive the governing balance equations by applying the principle of virtual work, which mandates that the sum of internal, external, and inertial forces acting on any material subdomain must vanish for any admissible virtual trial function. We employ Gurtin's microforce balance [19] to incorporate additional contributions of the entropy to the mechanical work as in the work of Forest and Amestoy [11]. Following the establishment of these balance equations, constitutive relations are derived from the first and second laws of thermodynamics, employing the well-established Coleman-Noll procedure [38]. Instead of working with a free energy, we work directly with the classical energy density. Then, for a rigid heat conductor, we obtain a fourth-order equation in terms of the entropy density that is capable of modeling scale-size effects for Fourier's law of heat conduction.

Solving the fourth-order differential equation that arises from the scale-size effects requires special attention in the choice of numerical methods. Standard finite element spaces based on $C^0$ continuity are insufficient due to the presence of higher-order derivatives in the governing equations. To ensure sufficient regularity in our approximation, we employ biquadratic B-splines [7], which naturally provide the required $C^1$ continuity across element boundaries. The tensor-product structure allows for a simple construction of multivariate basis functions for a finite element space, and for the computation of its derivatives in an univariate way. Therefore, to solve the governing equation and compute the effective thermal conductivity, we build a suitable bi-quadratic B-spline finite element approximation [20, 41].

The paper is organized as follows: In Section 2, we derive the model based on thermodynamic principles of continuum physics. In Section 3, we introduce a multiscale structure to the model and the associated internal lengths that capture the thermal scale-size effect. In Section 4, we deal with the numerical modeling of the aforementioned coefficients and challenges therein. Section 5 contains numerical simulations for different sizes of the representative volume elements of the volume fraction against the effective coefficients. Lastly, we reserve Section 6 for conclusions and perspectives on the work.

## 2. Thermodynamic modeling

We begin with some basic principles of continuum physics. Consider a body consisting of material points $X$ (a material point is identified with its position $X$ in the reference configuration [18]). The spatial position $x = \chi(X, t)$: here the function $\chi$, called the deformation, describes the motion of the body. Given a motion we distinguish between Grad (with upper-case G) and grad (with lower-case g) (likewise for Div and div): they refer to the spatial gradient operation (spatial divergence operation) with respect to $X$ and $x$ holding $t$ fixed (we will follow the same notation as Ogden [39].). Moreover, we use the superimposed dot to denote the material time derivative (i.e., with respect to $t$ holding $X$ fixed). We denote the deformation gradient
by $\mathbf{F} := \text{Grad}\boldsymbol{x}$, the velocity by $\boldsymbol{v} := \dot{\boldsymbol{x}}$, the mass density by $\rho$, the body forces by $\mathbf{b}$, the internal energy density by $e$, the internal entropy density by $s$, the absolute temperature by $T$, the heat flux by $\boldsymbol{q}$, and the heat source by $r$. Furthermore, we define $\mathcal{L}(\boldsymbol{v}) := \text{grad}\boldsymbol{v} = \dot{\mathbf{F}} \cdot \mathbf{F}^{-1}$ or when no confusion arises simply $\mathcal{L}$.

## 2.1. Conservation of mass

The conservation of mass has the following form,

$$ \dot{\varrho} + \varrho \operatorname{div} \boldsymbol{v} = 0. \qquad (2.1) $$

## 2.2. Balance of linear momentum

We will use the principle of virtual work (see, e.g., [15], [27]) to derive the general balance equations which asserts that the variation of the total elastic energy be equal to the virtual work exerted on the body by tractions, body forces, and acceleration forces acting on any subdomain $\Omega_t$ of the body. Define,

$$ W_{\text{int}}(\Omega_t; \boldsymbol{v}^*) := \int_{\Omega_t} \boldsymbol{\sigma} : \mathcal{L}(\boldsymbol{v}^*) \mathrm{d}\boldsymbol{v}, \qquad (2.2) $$

$$ W_{\text{ext}}(\Omega_t; \boldsymbol{v}^*) := \int_{\partial\Omega_t} \boldsymbol{t} \cdot \boldsymbol{v}^* \mathrm{d}\boldsymbol{a} + \int_{\Omega_t} \varrho \boldsymbol{b} \cdot \boldsymbol{v}^* \mathrm{d}\boldsymbol{v}, \qquad (2.3) $$

$$ W_{\text{acc}}(\Omega_t; \boldsymbol{v}^*) := - \int_{\Omega_t} \varrho \dot{\boldsymbol{v}} \cdot \boldsymbol{v}^* \mathrm{d}\boldsymbol{v}. \qquad (2.4) $$

The sum of the virtual work of internal, external, and acceleration forces acting on any subdomain $\Omega_t$ of the body for any virtual vector field $\boldsymbol{v}^*$ vanishes at each instant $t$ of the evolution:

$$ W_{\text{int}}(\Omega_t) = W_{\text{ext}}(\Omega_t) + W_{\text{acc}}(\Omega_t), \qquad (2.5) $$

for all $\Omega_t$ and for all admissible $\boldsymbol{v}^*$. Consequences of the virtual power principle are,

$$ \begin{aligned} \varrho \dot{\boldsymbol{v}} &= \operatorname{div}(\boldsymbol{\sigma}) + \varrho g && \text{in } \Omega_t, \\ \boldsymbol{t} &= \boldsymbol{\sigma} \cdot \boldsymbol{n} && \text{on } \partial\Omega_t. \end{aligned} \qquad (2.6) $$

## 2.3. Balance of angular momentum

We also need the balance of angular momentum,

$$ \frac{d}{dt} \int_{\Omega_t} \boldsymbol{x} \times \varrho \dot{\boldsymbol{x}} \mathrm{d}\boldsymbol{v} = \int_{\partial\Omega_t} \boldsymbol{x} \times [\boldsymbol{\sigma} \cdot \boldsymbol{n}] \mathrm{d}\boldsymbol{a} + \int_{\Omega_t} \boldsymbol{x} \times \varrho \boldsymbol{b} \mathrm{d}\boldsymbol{v}, \qquad (2.7) $$

which leads to the symmetry of $\boldsymbol{\sigma}$ when combined with the balance of linear momentum. In our case, the above argument is redundant since one of the consequences of one of the axioms of the method of virtual power is the symmetry of the stress (see the work of Maugin [27, Axiom 2.4, pg. 8–9]).

## 2.4. Entropy contribution to the mechanical work

As in the work of Forest and Amestoy [11], we indicate the micro-heterogeneous origin and scale size-dependence of the phenomena by the Greek letter $\chi$,

$$ W_{\text{int}}^{\chi}(\Omega_t; \dot{\boldsymbol{s}}^*) := \int_{\Omega_t} (a_s \dot{\boldsymbol{s}}^* + \boldsymbol{b}_s \cdot \text{grad} \, \dot{\boldsymbol{s}}^*) \mathrm{d}\boldsymbol{v}, \qquad (2.8) $$

$$ W_{\text{ext}}^{\chi}(\Omega_t; \dot{\boldsymbol{s}}^*) := \int_{\Omega_t} (a_s^p \dot{\boldsymbol{s}}^* + \boldsymbol{b}_s^p \cdot \text{grad} \, \dot{\boldsymbol{s}}^*) \mathrm{d}\boldsymbol{v} + \int_{\partial\Omega_t} a_s^S \dot{\boldsymbol{s}}^* \mathrm{d}\boldsymbol{a}, \qquad (2.9) $$
where $\dot{s}^*$ is the field of virtual rate of change of entropy density. The internal work of these defined quantities involves scalar and vector micro-stresses $a_s$, $b_s$ as well as prescribed external scalar and vector micro-forces $a_s^p$, $b_s^p$ as well as surface micro-traction $a_s^S$ (see [19], [11]).

Exploiting the principle of virtual work with respect to $\dot{s}^*$, we obtain the following independent (independent of the mechanical part) variational equations:

$$W_{\text{int}}^{\chi}(\Omega_t; s^*) = W_{\text{ext}}^{\chi}(\Omega_t; s^*), \quad \text{for all subdomains } \Omega_t \text{ and for all } \dot{s}^*, \qquad (2.10)$$

provided we have no micro-inertia effects for $s$. Hence, locally, we derive the following equations and associated boundary conditions,

$$\begin{aligned} \operatorname{div}(\mathbf{b}_s - \mathbf{b}_s^p) &= (a_s - a_s^p) && \text{in } \Omega_t \\ (\mathbf{b}_s - \mathbf{b}_s^p) \cdot \mathbf{n} &= a_s^S && \text{on } \partial\Omega_t \end{aligned} \qquad (2.11)$$

## 2.5. First law of thermodynamics

The balance law for the energy is then obtained,

$$\frac{d}{dt} \int_{\Omega_t} \rho e \, dv = W_{\text{ext}}(\Omega_t; \mathbf{v}) + W_{\text{acc}}(\Omega_t; \mathbf{v}) + W_{\text{ext}}^{\chi}(\Omega_t; \dot{s}) - \int_{\partial \Omega_t} \mathbf{q} \cdot \mathbf{n} \, da + \int_{\Omega_t} \rho r \, dv, \qquad (2.12)$$

or, equivalently, after taking the principle of virtual work into account,

$$\frac{d}{dt} \int_{\Omega_t} \rho e \, dv = W_{\text{int}}(\Omega_t, \mathbf{v}) + W_{\text{int}}^{\chi}(\Omega_t; \dot{s}) - \int_{\partial \Omega_t} \mathbf{q} \cdot \mathbf{n} \, da + \int_{\Omega_t} \rho r \, dv. \qquad (2.13)$$

The local equations then take the form,

$$\rho \dot{e} = \boldsymbol{\sigma} : \mathbf{L} + a_s \dot{s} + \mathbf{b}_s \cdot \operatorname{grad} \dot{s} - \operatorname{div}(\mathbf{q}) + \rho r. \qquad (\text{balance of energy})$$

## 2.6. Second law of thermodynamics

The second law of thermodynamics says that for any material subdomain $\Omega_t$,

$$\frac{d}{dt} \int_{\Omega_t} \rho s \, dv \geq - \int_{\partial \Omega_t} \frac{\mathbf{q}}{T} \cdot \mathbf{n} \, da + \int_{\Omega_t} \rho \frac{r}{T} \, dv. \qquad (\text{growth of entropy})$$

Combining the balance of energy with the growth of entropy inequality, we obtain the classical Clausius-Duhem inequality,

$$\rho(\dot{e} - T\dot{s}) \leq \boldsymbol{\sigma} : \mathbf{L} + a_s \dot{s} + \mathbf{b}_s \cdot \operatorname{grad} \dot{s} - \frac{1}{T} \mathbf{q} \cdot \operatorname{grad} T \qquad (\text{Clausius-Duhem ineq.})$$

For a smooth enough energy response function $e = \hat{e}(\mathbf{F}, s, \operatorname{grad} s)$ combined with the chain rule we can use the celebrated Coleman-Noll procedure (see [38]) and the relation $\dot{\mathbf{F}} \cdot \mathbf{F}^{-1}$ to obtain the following,

$$[\rho \frac{\partial \hat{e}}{\partial \mathbf{F}} - \boldsymbol{\sigma} \cdot \mathbf{F}^{-\top}] : \dot{\mathbf{F}} + [\rho \frac{\partial \hat{e}}{\partial s} - \rho T - a_s] \dot{s} + [\rho \frac{\partial \hat{e}}{\partial [\operatorname{grad} s]} - \mathbf{b}_s] \cdot \operatorname{grad} \dot{s} - \frac{1}{T} \mathbf{q} \cdot \operatorname{grad} T \leq 0, \qquad (2.14)$$

for all admissible processes<fnref n="1" />. From here, one can obtain the following constitutive laws,

$$\boldsymbol{\sigma} = \rho \frac{\partial \hat{e}}{\partial \mathbf{F}} \cdot \mathbf{F}^{\top}, \quad T = \frac{\partial \hat{e}}{\partial s} - \frac{1}{\rho} a_s, \quad \mathbf{b}_s = \rho \frac{\partial \hat{e}}{\partial [\operatorname{grad} s]} \qquad (2.15)$$

<footnote n="1">¹See the works of Noll and Coleman [38] or Coleman and Mizel [6] for a definition of what constitutes an *admissible process* based on the *Principle of Equipresence* defined in Truesdell and Toupin's work [46].</footnote>
Returning to the local equation of the balance of energy and using the chain rule for the energy density, we obtain,

$$0 = \left[ -\rho \frac{\partial \hat{e}}{\partial \mathbf{F}} + \boldsymbol{\sigma} \cdot \mathbf{F}^{-\top} \right] : \dot{\mathbf{F}} + \left[ -\rho \frac{\partial \hat{e}}{\partial s} + a_s \right] \dot{s} + \left[ -\rho \frac{\partial \hat{e}}{\partial [\mathrm{grad} \, s]} + \mathbf{b}_s \right] \cdot \mathrm{grad} \, \dot{s} - \mathrm{div}(\mathbf{q}) + \varrho r. \quad (2.16)$$

Using the constitutive laws from equation (2.15), we obtain the heat equation,

$$\varrho T \dot{s} = \mathrm{div}(\mathbf{q}) + \varrho r, \quad (2.17)$$

with the dissipation condition $-\frac{1}{T}\mathbf{q} \cdot \mathrm{grad} \, T \geq 0$.

## 2.7. Rigid heat conductor

Thus, for a rigid heat conductor<fnref n="2" /> the heat equation is,

$$0 = -\varrho T \dot{s} - \mathrm{div} \, \mathbf{q} + \varrho r \approx -\varrho T_0 \dot{s} - \mathrm{div} \, \mathbf{q} + \varrho r \quad \text{in } \Omega_t, \quad (2.18)$$

where $T_0$ is a strictly positive reference temperature. Furthermore, to satisfy the dissipation condition, we suppose that Fourier's law of heat conduction holds,

$$\mathbf{q} = -\mathbf{K} \cdot \nabla T \Rightarrow \mathbf{q} \cdot \nabla T \leq 0, \quad (2.19)$$

with $\mathbf{K}$ a second-order, positive definite heat conductivity tensor, which trivially guarantees that heat flows from hot to cold.

## 2.8. Approximation of the energy response function

The energy response function typically depends on the entropy and its gradient in a nonlinear way. However, if the energy response function is sufficiently smooth, we can approximate its behavior near equilibrium using a Taylor expansion. An equilibrium in this case means around the reference entropy $s_{\mathrm{ref}}$, which is the value of the entropy at absolute zero, which is usually unknown. In many problems, usually only the change or difference in entropy between different states is relevant, and, therefore, knowing the exact value of the entropy at absolute zero is not important. To further simplify this process, we introduce a scaled entropy variable, $\delta$, and define $s = \alpha \delta$, where $\alpha$ is a positive constant. By redefining the entropy in this scaled form, we perform a Taylor expansion up to order $\mathcal{O}(\alpha^3)$,

$$\begin{aligned} \hat{e}(s, \nabla s) = \hat{e}(s_{\mathrm{ref}}, \mathbf{0}) &+ \frac{\partial \hat{e}}{\partial s} (s_{\mathrm{ref}}, \mathbf{0}) (s - s_{\mathrm{ref}}) + \frac{\partial \hat{e}}{\partial [\nabla s]_i} (s_{\mathrm{ref}}, \mathbf{0}) [\nabla s]_i + \frac{1}{2} \frac{\partial^2 \hat{e}}{\partial s^2} (s_{\mathrm{ref}}, \mathbf{0}) (s - s_{\mathrm{ref}})^2 \\ &+ \frac{1}{2} \frac{\partial^2 \hat{e}}{\partial s \partial [\nabla s]_i} (s_{\mathrm{ref}}, \mathbf{0}) (s - s_{\mathrm{ref}}) [\nabla s]_i + \frac{1}{2} \sum_{i,j=1}^{3} \frac{\partial^2 \hat{e}}{\partial [\nabla s]_i \partial [\nabla s]_j} (s_{\mathrm{ref}}, \mathbf{0}) [\nabla s]_i [\nabla s]_j + \mathcal{O}(\alpha^3). \end{aligned} \quad (2.20)$$

We assume that we are working with a centro-symmetric material and, hence, any odd-order tensor is set to zero. Moreover, we set $\mathbf{b}_s^p = \mathbf{0}$ and $a_s^p = 0$ for expediting the presentation. Therefore,

$$\begin{aligned} \hat{e}(s, \nabla s) = \overbrace{\hat{e}(s_{\mathrm{ref}}, \mathbf{0})}^{\alpha} &+ \overbrace{\frac{\partial \hat{e}}{\partial s}(s_{\mathrm{ref}}, \mathbf{0})(s - s_{\mathrm{ref}})}^{\beta} + \frac{1}{2} \overbrace{\frac{\partial^2 \hat{e}}{\partial s^2}(s_{\mathrm{ref}}, \mathbf{0})(s - s_{\mathrm{ref}})^2}^{\gamma} \\ &+ \frac{1}{2} \sum_{i,j=1}^{3} \overbrace{\frac{\partial^2 \hat{e}}{\partial [\nabla s]_i \partial [\nabla s]_j} (s_{\mathrm{ref}}, \mathbf{0}) [\nabla s]_i [\nabla s]_j + \mathcal{O}(\alpha^3)}^{A_{ij}} \end{aligned} \quad (2.21)$$

<footnote n="2">²In the rigid heat conductor, there is no longer a difference between Grad and grad and thus we will replace it by $\nabla$. Moreover, the mass density $\varrho$ is constant.</footnote>
From here, one can compute,

$$ \frac{\partial \hat{e}(s, \nabla s)}{\partial s} = \frac{\partial \hat{e}}{\partial s} (s_0, \mathbf{0}) + \frac{\partial^2 \hat{e}}{\partial s^2} (s_0, \mathbf{0}) (s - s_0) + \mathcal{O}(\alpha^2), \qquad (2.22) $$

and

$$ \frac{\partial \hat{e}(s, \nabla s)}{\partial[\nabla s]} = \frac{\partial \hat{e}}{\partial[\nabla s]} (s_0, \mathbf{0}) + \frac{\partial^2 \hat{e}}{\partial[\nabla s]\partial[\nabla s]} (s_0, \mathbf{0}) : \nabla \nabla s + \mathcal{O}(\alpha^2). \qquad (2.23) $$

## 2.9. An additional length scale modeling thermal scale-size effects

We normalize the energy response function following a similar approach as in Nika and Muntean's work [37]. Let $\mathcal{A} := \max_{\forall t, x \in \Omega_t} |\mathbf{A}|$ and set $\mathbf{A}^* = \mathbf{A}/\mathcal{A}$. Moreover, we set $\mathcal{A} = \ell_s^2\gamma$. Then the scaled energy, call it $\hat{e}^*(s, \nabla s) = \hat{e}(s, \nabla s)/\gamma$.

Computing the relation between the temperature and entropy from the constitutive law for $T$, we obtain,

$$ T = \gamma(s - s_{ref}) - \ell_s^2 \sum_{i,j=1}^{3} \gamma A_{ij}^* \frac{\partial^2 s}{\partial x_i \partial x_j}. \qquad (2.24) $$

Using the above calculation for $T$ we can compute Fourier's law,

$$ q_i = - \sum_{j=1}^{3} K_{ij} \frac{\partial T}{\partial x_j} = - \sum_{i,j=1}^{3} \gamma K_{ij} \frac{\partial s}{\partial x_j} + \ell_s^2 \sum_{j,k,m=1}^{3} \gamma \overbrace{K_{ij} A_{km}^*}^{D_{ijkm}} \frac{\partial^3 s}{\partial x_j \partial x_k \partial x_m}. \qquad (2.25) $$

Therefore, the heat equation, for a homogeneous material, takes the form<fnref n="3" />,

$$ \varrho T_0 \frac{\partial s}{\partial t} = \sum_{i,j=1}^{3} \gamma K_{ij} \frac{\partial^2 s}{\partial x_i \partial x_j} + \ell_s^2 \sum_{i,j,k,m=1}^{3} \gamma K_{ij} A_{km} \frac{\partial^4 s}{\partial x_i \partial x_j \partial x_k \partial x_m} + \varrho r. \qquad (2.26) $$

We remark that this non-classical thermodynamic derivation of the heat equation in terms of the entropy density leads to equation (2.26), which is a fourth-order linear parabolic partial differential equation. A similar structured evolutionary equation was derived by Nika [33] within the GENERIC framework for the absolute temperature and using a classical gradient flow through an Onsager system with quadratic mobility. There, the main idea for the derivation (see also [36]) was to postulate that the heat flux depends linearly on the temperature gradient but also on high-order temperature derivatives (see [6, 40, 5]). The reader is encouraged to consult the work of Matthes et al. [26] for the use of different gradient flow metrics with second-order derivatives for a derivation of similar structured equations with linear mobility.

# 3. Multiscale analysis of a prototype diffusion with scale-size thermal effects

## 3.1. The multiscale structure of the heterogeneous medium

The study of heat flow evolution in composite materials holds significant importance for applied scientists, as the ability to regulate constituent proportions and their absolute dimensions enables precise control over thermal properties to achieve targeted outcomes. For example, strategic material design can facilitate spatially varying diffusion rates—slower in specific regions and faster in others—thereby passively directing heat flow without external intervention. In such systems, scale-dependent effects critically influence performance, as evidenced by recent experimental investigations of thermal behavior in heterogeneous media [9, 10]. To model these phenomena, we employ the governing equation (2.26) within a periodic framework.

<footnote n="3">³We remark that when $\mathbf{K} = \kappa \mathbf{I}$ and $\mathbf{A} = \alpha \mathbf{I}$, with $\kappa > 0$ and $\alpha > 0$, then 2.26 takes the form $\varrho T_0 \frac{\partial s}{\partial t} = \kappa \Delta s - \alpha \Delta^2 s$ which is the classical structure of a Cahn-Hilliard equation, as pointed out in [11].</footnote>
Specifically, the composite material occupies a domain $\bar{\Omega} \subset \mathbb{R}^3$, characterized by periodicity $\ell$, macroscopic length L, and a uniformly Lipschitz exterior boundary $\partial\Omega$ with outward unit normal $n$. The microstructure is represented by a periodic unit cell Y, partitioned into a simply connected subdomain and its complement, corresponding to distinct material phases (see Fig. 2). The tensors $K$ and $D$ are assumed to be

<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic of the periodic microstructure of the composite in three dimensions.</figcaption></figure>

$Y := (-\ell/2, \ell/2]^3$ periodic, with the periodic cell Y characterizing the body $\Omega$, and piece-wise constant in each phase. They stand for the thermal conductivity and spatial retardation, respectively. We consider a particular form of equation (2.26),

$$ \varrho T_0 \frac{\partial s}{\partial t} = \operatorname{div} (\gamma K \cdot \nabla s - \operatorname{div} (\gamma \ell_s^2 D : \nabla \nabla s)) + \varrho r \quad \text{in } (0, T) \times \Omega. \quad \text{(enhanced heat equation)} $$

## 3.2. Scaling of the variables

We normalize all variables and material properties using the following transformations:

$$ t^* = \frac{t}{\tau}, \quad x^* = \frac{x}{L}, \quad s^*(x^*) = \frac{s(x)}{L}, \quad \mathcal{K}K^* = K \quad (3.1) $$

where L denotes the macroscopic characteristic length of the domain and $\mathcal{K} := \max_{z \in Y} |K(z)|$ the characteristic thermal conductivity value associated with tensor $K$. Two distinct time scales emerge from this normalization:

$$ \tau = L^2/\mathcal{K} \quad (\text{macroscopic}) \quad \text{and} \quad \tau = \ell^2/\mathcal{K} \quad (\text{microscopic}), \quad (3.2) $$

corresponding to diffusion timescales at their respective spatial scales. As our focus lies in deriving effective macroscopic equations, we consider only the slower macroscopic timescale (see [28]).

Applying these normalizations to (2.26) yields:

$$ \varrho T_0 \gamma^{-1} \frac{\partial s^*}{\partial t^*} = \operatorname{div}^* \left( K^* \cdot \nabla^* s^* - \left(\frac{\ell_s}{L}\right)^2 \operatorname{div}^* (D^* : \nabla^* \nabla^* s^*) \right) + \varrho \gamma^{-1} r \quad \text{in } (0, T) \times \Omega, \quad (3.3) $$

where the scaled coefficients $K^*$ and $D^*$ are $Y^*$ periodic where $Y^* := (\ell/L)Y$ with $Y := (-1/2, 1/2]^3$. For notational simplicity, we subsequently omit the asterisk (*) superscripts.

Assuming scale separation between the microscopic ($\ell$) and macroscopic (L) lengths through the ratio $\varepsilon := \ell/L \ll 1$, we consider the homogenization limit as $\varepsilon \to 0$. While classical homogenization operates with these two fundamental scales, the formulation used within this article introduces an additional non-classical consideration: the relative magnitude of the slip length $\ell_s$ compared to either $\ell$ or L. This distinction becomes crucial in multiscale modeling of periodic heterogeneous rigid heat conductors beyond conventional approaches. It is important to point out that the term slip length arises from the theory of rarefied gas flow. For a rarefied gas at the slip regime, namely when the Knudsen number is $0.001 < K_n < 0.3$. The
prototypical experiment conducted is to consider a steady heat flux through argon gas between two infinite parallel plates, which are set at unit distance in the slip-regime, as defined above. In such an experiment, the classical Fourier's law of heat conduction cannot correctly describe the temperature field in a rarefied gas at the slip regime. The above failure to correctly describe the temperature field in the slip regime has led to the speculation that heat flux in a rarefied gas may additionally depend on higher-order temperature derivatives. Thus, it is possible to obtain higher-order corrections to the flux in a Taylor series in terms of the Knudsen number. When the Knudsen number approaches zero, the higher-order derivative terms can be neglected, and Fourier's law is recovered. However, at the slip regime, if the heat flux varies quadratically with the Knudsen number, then the contribution of the higher-order terms may not be a negligibly small quantity. The length $l_s$ represents the quadratic variation of the heat flux with respect to the Knudsen number.

### 3.3. Slip length comparable to the size of the heterogeneity

When the slip length is comparable to the size of the heterogeneity, $l_s/l \sim 1$, then this implies that,

$$l_s/L \sim \varepsilon. \qquad (\text{scaling})$$

Hence, with $\varepsilon \ll 1$ we remain in the classical continuum regime since as $\varepsilon$ goes to zero this implies that the corresponding Knudsen number is below $10^{-3}$ and therefore the heat flow can be treated as a classical continuum. Thus, the enhanced heat equation under this scaling becomes

$$\varrho T_0 \gamma^{-1} \frac{\partial s_{\varepsilon}}{\partial t} = \text{div} \left( \mathbf{K} \left( \frac{\mathbf{x}}{\varepsilon} \right) \cdot \nabla s_{\varepsilon} - \varepsilon^2 \text{div} \left( \mathbb{D} \left( \frac{\mathbf{x}}{\varepsilon} \right) : \nabla \nabla s_{\varepsilon} \right) \right) + \varrho \gamma^{-1} r \quad \text{in } (0, T) \times \Omega, \quad (3.4)$$

The microscopic problem is then characterized by the following equation and boundary & initial conditions,

$$\varrho T_0 \gamma^{-1} \frac{\partial s_{\varepsilon}}{\partial t} = \text{div} \left( \mathbf{K} \left( \frac{\mathbf{x}}{\varepsilon} \right) \cdot \nabla s_{\varepsilon} - \varepsilon^2 \text{div} \left( \mathbb{D} \left( \frac{\mathbf{x}}{\varepsilon} \right) : \nabla \nabla s_{\varepsilon} \right) \right) + \varrho \gamma^{-1} r \quad \text{in } (0, T) \times \Omega, \quad (3.5)$$

$$[\varepsilon^2 \mathbb{D}(\frac{\mathbf{x}}{\varepsilon}):\nabla\nabla s_{\varepsilon}] : [\mathbf{n} \otimes \mathbf{n}] = 0 \qquad \text{on } (0, T) \times \partial\Omega, \quad (3.6)$$

$$s_{\varepsilon}(t, \mathbf{x}) = 0 \qquad \text{on } (0, T) \times \partial\Omega, \quad (3.7)$$

$$s_{\varepsilon}(0, \mathbf{x}) = s_{\varepsilon}^{0}(\mathbf{x}) \qquad \text{in } \Omega. \quad (3.8)$$

where $r$ is a given source that belongs in $L^2(\Omega)$. We have prescribed a homogeneous Dirichlet boundary condition in (3.7), however, this is no longer sufficient. We require, additionally, to prescribe the normal double heat flux in condition (3.6). It plays, in part, the role of a natural boundary condition. Its dual consists of fixing $\nabla s_{\varepsilon} \cdot \mathbf{n}$ on the boundary. In direct analogy with second-gradient elasticity, such a term is referred to as normal double-traction (see, in particular, [15], [16], [31]).

### 3.4. Homogenization via two-scale expansion for the steady state enhanced heat equation

This section aims to derive effective equations to link and describe the microscopic and the macroscopic behavior of the high-order system at steady-state via two-scale expansion [3, 42, 1, 28]. Let us assume the tensors $\mathbf{K}$ and $\mathbb{D}$ to be bounded and uniformly elliptic to ensure existence and uniqueness of the solution to the high-order system [37]. The effective properties also depend on the absolute size of the heterogeneity [12]. This is made clear by the chosen scaling, indicating that the microscopic thermal effects and the size of the heterogeneity have comparable magnitudes, ensuring full scale separation$^4$.

<footnote>$^4$A different scaling choice could lead to the various impacts of the higher gradient in terms of both relevance and length scale [34], [35]. By choosing, for example, $l_s/L \sim 1$, one would attain partial scale separation between effective properties and the heterogeneity size. Such intrinsic lengths related to scale-dependent thermal effects have been observed in experiments, see, e.g., [9, 10]</footnote>
Consider a domain $\Omega \in \mathbb{R}^n$ with periodic microstructure $Y$ of period $\varepsilon$. Then, in the steady-state case, equations (3.5) to (3.8) become

$$ 
\begin{aligned}
& -\operatorname{div}\left(\mathbf{K}\left(\frac{\mathbf{x}}{\varepsilon}\right) \cdot \nabla s_{\varepsilon}-\varepsilon^{2} \operatorname{div}\left(\mathbb{D}\left(\frac{\mathbf{x}}{\varepsilon}\right) \cdot \nabla \nabla s_{\varepsilon}\right)\right)=\varrho \gamma^{-1} r \quad \text { in } \Omega, \\
& \left[\varepsilon^{2} \mathbb{D}\left(\frac{\mathbf{x}}{\varepsilon}\right) \cdot \nabla \nabla s_{\varepsilon}\right]:\left[\mathbf{n} \otimes \mathbf{n}\right]=0 \quad \text { on } \partial \Omega, \\
& s_{\varepsilon}(\mathbf{x})=0 \quad \text { on } \partial \Omega,
\end{aligned}
\quad (3.9) $$

where $\mathbf{K}$ and $\mathbb{D}$ are bounded and uniformly elliptic tensors. To homogenize the above system, we assume sufficient regularity and scale separation on all of the involved functions, which can be expressed in terms of the slow variable $\mathbf{x}$ and the fast variable $\mathbf{y} = \mathbf{x}/\varepsilon$.

In order to obtain the homogenized equation, we two-scale expand the solution, $s_\varepsilon$, as follows:

$$ s_{\varepsilon}(\mathbf{x}) = \sum_{i=0}^{+\infty} \varepsilon^{i} s^{(i)}(\mathbf{x}, \mathbf{y}), \quad \text{with} \quad \mathbf{y} = \mathbf{x}/\varepsilon. \quad (3.10) $$

Moreover, we also express the gradient as $\nabla \rightsquigarrow \nabla_x + \varepsilon^{-1}\nabla_y$ and the second gradient as $\nabla \nabla \rightsquigarrow \nabla_x \nabla_x + 2\varepsilon^{-1}\nabla_x \nabla_y + \varepsilon^{-2}\nabla_y \nabla_y$ to include the effects of both scales. Using the above two-scale expansion, the derivatives of the entropy are

$$ \nabla s_{\varepsilon} = \varepsilon^{-1} \nabla_y s^{(0)} + \varepsilon^0 (\nabla_x s^{(0)} + \nabla_y s^{(1)}) + \varepsilon^1 (\nabla_x s^{(1)} + \nabla_y s^{(2)}) + \dots, \quad (3.11) $$

and

$$ 
\begin{aligned}
\nabla \nabla s_{\varepsilon} = & \varepsilon^{-2} \nabla_y \nabla_y s^{(0)} + \varepsilon^{-1} (2 \nabla_x \nabla_y s^{(0)} + \nabla_y \nabla_y s^{(1)}) \\
& + \varepsilon^0 (\nabla_x \nabla_x s^{(0)} + 2 \nabla_y \nabla_x s^{(1)} + \nabla_y \nabla_y s^{(2)}) + \dots
\end{aligned}
\quad (3.12) $$

For the flux, the two-scale formula translates as

$$ \operatorname{div} \boldsymbol{q}_{\varepsilon}=\varepsilon^{0} \operatorname{div}_{y} \boldsymbol{q}^{0}+\varepsilon^{-1}\left(\operatorname{div}_{x} \boldsymbol{q}^{0}+\operatorname{div}_{y} \boldsymbol{q}^{1}\right)+\varepsilon^{0}\left(\operatorname{div}_{x} \boldsymbol{q}^{1}+\operatorname{div}_{y} \boldsymbol{q}^{2}\right)+\dots \quad (3.13) $$

We plug the two-scale expansion in System (3.9) and study the system's behavior at different length scales.

Matching the orders for $\varepsilon^{-2}$, we attain

$$ 
\begin{aligned}
& -\operatorname{div}_y \left( \mathbf{K}(\mathbf{y}) \nabla_y s^{(0)} - \operatorname{div}_y (\mathbb{D}(\mathbf{y}):\nabla_y \nabla_y s^{(0)}) \right) = 0 && \text{in } Y, \\
& s^{(0)}(\mathbf{x}, \mathbf{y}) && \text{periodic in } Y.
\end{aligned}
\quad (3.14) $$

Given the assumptions on the coefficients, the above system has a unique solution in the fast variable $\mathbf{y}$ by the Lax-Milgram lemma<fnref n="5" />. Given that any function of $\mathbf{x}$ is also a solution to the above system, we have $s^{(0)}(\mathbf{x}, \mathbf{y}) = s^{(0)}(\mathbf{x})$.

For $\varepsilon^{-1}$, we have the microscopic problem,

$$ 
\begin{aligned}
& -\operatorname{div}_y \left( \mathbf{K}(\mathbf{y}) (\nabla_x s^{(0)} + \nabla_y s^{(1)}) - \operatorname{div}_y (\mathbb{D}(\mathbf{y}):\nabla_y \nabla_y s^{(1)}) \right) = 0 && \text{in } Y, \\
& s^{(1)}(\mathbf{x}, \mathbf{y}) && \text{periodic in } Y.
\end{aligned}
\quad (3.15) $$

From the above structure of the equation, we can see that the unknown $s^{(1)}$ depends linearly on $\nabla_x s^{(0)}$. By virtue of this linearity we can decompose the contribution of $s^{(1)}$ as follows,

$$ s^{(1)}(\mathbf{x}, \mathbf{y}) = \sum_{i=1}^{3} w^{i}(\mathbf{y}) \frac{\partial s^{(0)}}{\partial x_{i}}(\mathbf{x}) + c(\mathbf{x}), \quad (3.16) $$

<footnote n="5">⁵One has to work on the Sobolev space $H^2(Y)$ with $Y$ periodic functions.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3: Four possible microscopic periodic unit cells $Y_m$, where $m^2$ indicates the number of microstructures in the cell, which is related to the “absolute amount of constituents” in each cell.</figcaption></figure>

where $w^i(\mathbf{y})$ are called corrector terms and $c(\mathbf{x})$ is an arbitrary additive constant (function of $\mathbf{x}$). Rewriting system (3.15) in terms of the corrector terms yields the cell problems. For $i = 1, 2, 3$ and $e_i \in \mathbb{R}^3$ unit vector in $i^{th}$ direction, the $i^{th}$ cell problem is,

$$ -\operatorname{div}_y (\mathbf{K}(\mathbf{y})(e_i + w^i(\mathbf{y})) - \operatorname{div}_y (\mathbb{D}(\mathbf{y}):\nabla_y \nabla_y w^i(\mathbf{y}))) = 0 \qquad \text{in } Y, \tag{3.17} \\ w^i(\mathbf{y}) \qquad \text{periodic in } Y. $$

Note that the corrector terms $w^i$ solving the respective cell problems at the microscopic level contain relevant information on the behavior of the heterogeneous medium. Moreover, a direct consequence of the scaling choice $\ell_s/\ell \sim 1$ is the higher-order term appearing only in the microscopic equation.

Matching the orders for $\varepsilon^0$, the macroscopic problem is,

$$ -\operatorname{div}_x (\mathbf{K}^{\text{eff}} \nabla_x s^{(0)}) = \varrho \gamma^{-1} r \qquad \text{in } \Omega, \tag{3.18} \\ s^{(0)}(\mathbf{x}) = 0 \qquad \text{on } \partial\Omega, $$

where $\mathbf{K}^{\text{eff}}(\mathbf{y})$ is the tensor describing the effective conductivity of the homogenized medium in relation to the cell problem. Namely, for $i, j = 1, 2, 3$, the $ij^{th}$ entry of the effective conductivity tensor $K_{ij}^{\text{eff}}$ is defined as,

$$ K_{ij}^{\text{eff}} = \sum_{k=i}^{3} \frac{1}{|Y|} \int_{Y} K_{ik}(\mathbf{y}) \left( \delta_{ik} + \frac{\partial w^j}{\partial y_k} \right) d\mathbf{y}. \tag{3.19} $$

The information on microscale behaviors is transferred to the macroscopic equation through the effective conductivity tensor; thus, the macroscopic equation does not depend explicitly on the higher-order derivative. We stress at this point that although the formula for the effective conductivity in our scaling regime, assumed in our scaling, looks like the classical one, it is not. The difference is in how the corrector solutions $w^j$ are computed, namely, one must account for higher-gradient terms appearing in the cell problem (3.17).

## 4. Numerical modeling

To numerically compute the effective thermal conductivity $\mathbf{K}^{\text{eff}}$ in a two-dimensional setup, we need to solve the cell problems (3.17), taking into account the higher-order gradient and the consequent regularity required of the test space. Figure 3 shows four possible microscopic periodic unit cells $Y_m$, each containing $m^2$ microstructures, which are scaled versions of the base unit cell $Y_1$. Experimental results in a similar context [32] show that the second-gradient term affects the effective properties depending on the absolute size of the microstructure, that is, in the periodic cell $Y_1$, the second-gradient effects are dominant ($\ell_s/\ell \sim 1$), while for large values of $m$ they are negligible ($\ell_s/\ell \sim \varepsilon$). Experimentally, it has been determined that a
<figure><img src="image_5.png" /><figcaption>Figure 4: Support of two B-splines on the domain $[0, 1]^2$ uniformly divided into 8 × 8 square elements of side $h = 1/8$. The squares highlighted in cyan in the central part of the domain indicate the support of $B_{l,h}$, while the ones highlighted in pink indicate the support of $B_{k,h}$ across the periodic boundary.</figcaption></figure>

ratio of 1/56 between micro- and macroscales ($\ell_s/\ell \sim 1/56$) is sufficient to achieve strict scale separation and ignore second-gradient effects [32].

Here, we focus on the case $K = \kappa I$, where $\kappa$ is a constant. Moreover, for any twice differentiable function $\phi$, we relate the $ijkl$th entry of the spatial retardation term $D$ to the conductivity $K$ as

$$D_{ijkl} \partial_{y_k} \partial_{y_l} \phi = \gamma_m^2 K_{ik} \partial_{y_k} \partial_{y_l} \phi = \gamma_m^2 \kappa \delta_{ik} \partial_{y_k} \partial_{y_l} \phi = \gamma_m^2 \kappa \partial_{y_i} \partial_{y_j} \phi, \quad (4.1)$$

for $i, j, k, l = 1, 2$, where $I$ is the identity tensor, and $\gamma_m$ is a length-scale parameter depending on the number of microstructures in the periodic unit cell. The length-scale parameter is such that, for small values of $\gamma_m$, the high-order term becomes negligible and the cell problem reduces to the classical elliptic problem, that is, we recover the classical homogenization theory with scale separation between the effective properties and the heterogeneity size ($\ell_s/\ell \sim \epsilon$).

The variational problems of cell problems reads: find $w^i$ for $i = 1, 2$ and for every appropriate test function $\phi$, such that

$$\int_{Y_m} \kappa \nabla w^i \cdot \nabla \phi + \int_{Y_m} \gamma_m^2 \kappa \nabla \nabla w^i : \nabla \nabla \phi = - \int_{Y_m} \kappa e_i \cdot \nabla \phi, \quad (4.2)$$

and $w^i$ are periodic in $Y_m$. In the remainder of this section, we construct a suitable B-spline finite element method [20, 41] to solve cell problems (3.17) and compute the effective thermal conductivity $K^{\text{eff}}$.

The cell problem necessitates a numerical discretization scheme with sufficient smoothness across element boundaries. To meet this requirement, we employ bi-quadratic B-splines, which provide $C^1$ continuity across elements and are well suited for solving fourth-order PDEs on structured grids. These B-splines offer several advantages: they possess compact support and facilitate efficient implementation, and the tensor-product structure allows for straightforward construction of multivariate basis functions from univariate components.

We discretize the computational domain into a uniform grid, and the spline basis functions are centered on the element centers. Special care is taken to ensure consistency across the periodic boundaries, where periodicity is enforced by identifying degrees of freedom across opposite boundaries. That is, each bi-quadratic B-spline basic function spans a 3 × 3 patch of elements; Figure 4 illustrates how the periodic nature of the domain $[0, 1]^2$ affects the support of the B-splines on a periodic domain discretized into 8 × 8 elements. While this tensor product structure does not fully explore the potential of multivariate B-spline features [7], in our case, it provides a simple and effective way to construct basis functions and compute
their derivatives in a univariate manner. To build the finite element basis, we consider only a part $Y^{(m)}$ of the unit cell that contains only a single microstructure, that is, $Y^{(m)} = Y_1$ if $m = 1$, or a portion of $Y_m$—a scaled version of $Y_1$—if $m > 1$. Note that the side length of $Y^{(m)}$ is $1/m$, and the size effects enter as the different terms in the variational problem scale differently. To implement a B-spline based finite element solver for equation (4.2), we divide the domain $Y^{(m)}$ into $N = n \times n$ square elements of side $h = 1/(nm)$. To keep notation simple, we let $\varphi_1^{(h)}, \varphi_2^{(h)}, \dots, \varphi_N^{(h)}$ denote the basis functions on $Y^{(m)}$, the $1/m$ scaled version of $Y_1$.

We employ a classical Galerkin approximation where we set both the test and trial function spaces as $V_h = \text{span}\{\varphi_1^{(h)}, \varphi_2^{(h)}, \dots, \varphi_N^{(h)}\}$. That is, the solutions $w_h^i$, $i = 1, 2$, are approximated by,

$$w_h^i = \sum_{k=1}^{N} d_k^i \varphi_k^{(h)}, \qquad (4.3)$$

where $d_k^i$, $i = 1, 2$, $k = 1, \dots, N$ are the coefficients to be determined. Unlike standard Lagrangian basis functions, the coefficients $d_k^i$ do not correspond to the nodal values of $w_h^i$—however, they can be recovered efficiently [14]. For simplicity, we approximate, the conductivity $\kappa$ by an element-wise constant function $\kappa_h$ in the numerical experiments below.

Let $d^i \in \mathbb{R}^N$, $i = 1, 2$, be the column vectors holding the unknown coefficients and let $1 = (1, 1, \dots, 1)^T \in \mathbb{R}^N$, then the algebraic form variational problem (4.2) is,

$$\begin{bmatrix} S + D & \mathbf{1} \\ \mathbf{1}^T & 0 \end{bmatrix} \begin{bmatrix} d^i \\ \lambda^i \end{bmatrix} = \begin{bmatrix} f^i \\ 0 \end{bmatrix}, \qquad (4.4)$$

where $\lambda^i$ is a Lagrange multiplier corresponding to the periodic boundary condition, and the stiffness matrix $S$, the spatial retardation matrix $D$, and the right-hand side vector $f^i$ have entries,

$$s_{kl} = \int_Y \kappa_h \nabla \varphi_k^{(h)} \cdot \nabla \varphi_l^{(h)}, \quad d_{kl} = \int_Y \gamma_m^2 \kappa_h \nabla \varphi_k^{(h)} : \nabla \varphi_l^{(h)}, \qquad (4.5)$$

and

$$f_k^i = - \int_Y \kappa_h e_i \cdot \nabla \varphi_k^h, \quad i = 1, 2, \qquad (4.6)$$

respectively. It is worth noting that both cell problems share the same system matrix, which is real and symmetric. This allows for a single matrix factorization to be reused across both solves, as well as in potential future adjoint computations. Such reuse significantly reduces computational cost and is essential for future large-scale simulation and optimization runs.

## 5. Numerical simulations

We numerically investigate the effects of the spatial retardation term on the effective thermal conductivity for various microstructure shapes and sizes. The length-scale parameter $\gamma_m = 1/m^2$ characterizes the absolute size of the constituents in the unit cell $Y_m$. All simulations are performed on a square domain $Y^{(m)} = [0, 1]^2$, which is discretized into a uniform grid with $n = 512 \times 512$ elements.

Figure 5 displays the microstructures used in the numerical experiments. For every microstructure, the yellow region indicates a material with high thermal conductivity, characterized by the parameter $\kappa_c = 250$, while the blue region indicates a material with low thermal conductivity, characterized by the parameter $\kappa_i = 0.2$. Figure 6 shows, from left to right, the different domains $Y_m$, $m = 1, \dots, 6$, each incorporating the circular microstructure shown in Figure 5. Note that, for $m \to \infty$, the problem reduces to the classical homogenization scenario, and we denote the domain by $Y^{(\infty)}$.

We solve the cell problems (4.2) on both the circular and the square microstructures, highlighting the effects of different constituent sizes. Figure 7 shows the solutions $w^1$ and $w^2$ of the cell problems on the
<figure><img src="image_6.png" /><figcaption>Figure 5: Microstructures used for the numerical experiments. The yellow region represents the high-conductivity material ($\kappa_c = 250$), while the blue region is low-conductivity material ($\kappa_i = 0.2$).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Periodic unit cells $Y_m$, for $m = 1, \dots, 6$, sharing the same microstructural geometry but differing in characteristic size; each unit cell contains $m \times m$ scaled versions of the underlying microstructure. The yellow portion indicates a high conductivity denoted by $\kappa_c$, while the blue portion indicates a low conductivity denoted by $\kappa_i$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Local solutions $w^i$, $i = 1, 2$ computed from equation (4.2) for the circular microstructure and shown on microstructure cells $Y^{(m)}$ for $m = 1, \dots, 6$ illustrating the effect of the absolute size of the constituents on the coefficients.</figcaption></figure>

microstructure cells $Y^{(m)}$ for the circular microstructure, where the number of microstructures per side is $m = 1, \dots, 6$. Analogously, Figure 8 shows the solutions $w^1$ and $w^2$ of the cell problems on the microstructure cell $Y^{(m)}$, for the square microstructure and $m = 1, \dots, 6$. In both cases, the solutions are shown on the microstructure cells $Y^{(m)}$ to highlight the effect of different absolute sizes of constituents on the solutions.

To examine the behavior of the effective properties as the absolute size of constituents decreases, we compare the effective conductivity of the scale-dependent problems model to that computed using classical homogenization theory. Figure 9 shows the relative difference in $K_{11}^{\text{eff}}$ between the scale-dependent case ($m = 1, \dots, 100$), and the classical case ($m \to \infty$), for both the circular and the square microstructures. As expected, the relative difference in $K_{11}^{\text{eff}}$ decreases as the number of microstructures increases, indicating convergence towards the classical case as the absolute sizes of constituents decrease. In contrast, for both the circular and square microstructure, $K_{12}^{\text{eff}}$ appears to be independent of $m$. More precisely, for a fixed microstructure geometry, the relative difference between $K_{12}^{\text{eff}}$ computed on microstructure cells $Y^{(m)}$ and the classical homogenized case is smaller than the numerical error introduced by the finite element solver at
<figure><img src="image_9.png" /><figcaption>Figure 8: Local solutions $w^i$, $i = 1, 2$ computed from equation (4.2) for the square microstructure and shown on microstructure cells $Y^{(m)}$ for $m = 1, \dots, 6$ characterizing the effect of the absolute size of the constituents on the coefficients.</figcaption></figure>

the current resolution (512 $\times$ 512 elements).

We study the effect of microstructure size relative to the absolute size of the constituents, by simulating microstructure cells $Y^{(m)}$, for selected values of $m$, with varying volume fractions. Recall that for the extreme cases of 0% and 100% volume fraction, the problem reduces to a uniform conductivity of $\kappa_i$ and $\kappa_c$, respectively. Figure 10 and Figure 11 show how the solution $w^1$ varies with volume fraction for circular and square microstructures, respectively. More precisely, the figures show the solution $w^1$ computed on microstructure cells $Y^{(m)}$, for $m = 1, \dots, 5$, and for the classical homogenization case ($m \to \infty$) at volume fractions of 5%, 40%, and 75%.

Figure 12 shows the effective coefficients $K_{11}^{\text{eff}}$ plotted across different volume fractions (from 5% to 75%) for circular and square microstructures. Each plot contains the effective coefficients computed on the microstructure cells $Y^{(m)}$, for $m = 1, \dots, 5$, and for the classical homogenization case. The behavior of $K_{11}^{\text{eff}}$ for different volume fractions converges towards the coefficients computed in the classical homogenization case as we increase $m$. This aligns with intuition, as $m \to \infty$, the microstructure becomes vanishingly small,

<figure><img src="image_10.png" /><figcaption>Figure 9: Relative difference in $K_{11}^{\text{eff}}$ between scale-dependent case (with $m = 1, \dots, 100$ microstructures per side) and the classical homogenization ($m \to \infty$), for both circular and the square microstructure configurations.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Local solutions $w^1$ computed for the circular microstructure with varying volume fractions (5%, 40%, and 75%) on microstructure cells $Y^{(m)}$, $m = 1, \dots, 5$ and for the classical homogenization case ($m \to \infty$).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Local solutions $w^1$ computed for the square microstructure with varying volume fractions (5%, 40%, and 75%) on microstructure cells $Y^{(m)}$, $m = 1, \dots, 5$ and for the classical homogenization case ($m \to \infty$).</figcaption></figure>

and classical homogenization applies. In contrast, $K_{11}^{\text{eff}}$, depending on scale-size effects, $K_{12}^{\text{eff}}$ appears to be independent of the absolute size of the constituents. This behavior might be a consequence of the modeling choice in (4.1), with $K = \kappa I$.
<figure><img src="image_13.png" /><figcaption>Figure 12: Effective coefficients $K_{11}^{\text{eff}}$ and $K_{12}^{\text{eff}}$ plotted across different volume fractions for microstructure cells with varying length-scales and geometric configurations (circular and square), respectively. For both geometries, scale-effects influence $K_{11}^{\text{eff}}$ but not $K_{12}^{\text{eff}}$.</figcaption></figure>

# 6. Concluding remarks

We homogenized an enhanced heat equation, taking into account thermal scale-size effects, provided a numerical framework for the simulation of the derived homogenized effective coefficients based on quadratic splines in order to account for the higher gradients, and computed several examples for different-sized representative volume elements of the volume fraction versus the effective coefficients. The framework's usability is not limited to analyzing scale-dependent behavior; it also enables the use of material distribution topology optimization. By leveraging the effective coefficients derived here, one can design microgeometries with desired thermal properties, such as anisotropic conductivity or localized heat control. While localized heat control would require spatially varying microstructures across the macroscale—potentially drawn from a continuous catalog of microstructure designs, such as the circular and square inclusions for varying volume fractions—the current framework provides the tools to handle such configurations. This opens up promising
opportunities for designing thermal metamaterials with engineered heat transport characteristics.

A similar structured equation was derived by Nika [33] within the GENERIC framework for the absolute temperature and using a classical gradient flow through an Onsager system with quadratic mobility. The reader is also encouraged to consult the work by Matthes et al. [26] for the use of different gradient flow metrics with second-order derivatives for a derivation of similar structured equations with linear mobility. The main idea for the derivation by Nika [33] (see also Nika and Muntean [36]) was to postulate that the heat flux depends linearly on the temperature gradient but also on high-order temperature derivatives [5, 6, 40]. Pan et al. [40] argued—albeit in the context of temperature—that such a higher-gradient term reflects the influence of the Knudsen layer in rarefied gases in the *slip-regime*. Hence, Fourier's law of heat conduction was not adhered to, and the dissipation due to heat needed to be imposed to have a thermodynamically consistent model. In the current approach, the above scenario is evaded by working directly with the entropy density and Fourier's law of heat conduction; hence, the dissipation inequality due to heat was automatically satisfied.

## Code availability

The code is available from the corresponding author upon reasonable request.

## Replication of results

All the details necessary to reproduce the results have been defined in the paper.

## Funding

GN gratefully acknowledges the financial support from the Swedish Research Council (VR) (Grant no. 2024-04069), the Knowledge Foundation (project nr. KK 2020-0152). EW gratefully acknowledges the financial support from the Swedish Research Council (VR) (Grant no. 2022-03783).

## Acknowledgements

GN thanks the group of Michael Stingl at FAU Erlangen-Nürnberg for their hospitality where part of this work was completed.

## Author contributions

Conceptualization: Mario Setta, Eddie Wadbro, Grigor Nika; Funding acquisition: Eddie Wadbro, Grigor Nika; Methodology: Mario Setta, Eddie Wadbro, Grigor Nika; Software: Mario Setta, Eddie Wadbro; Visualization: Mario Setta, Eddie Wadbro, Grigor Nika; Writing – original draft: Mario Setta, Eddie Wadbro, Grigor Nika; Writing – review and editing: Mario Setta, Eddie Wadbro, Grigor Nika;

## References

[1] N. Bakhvalov and G. Panasenko. *Homogenisation: averaging processes in periodic media: mathematical problems in the mechanics of composite materials*. Kluwer Academic Publishers, 1989.

[2] J. Bauer, L. R. Meza, T. Schaedler, R. Schwaiger, X. Zheng, and L. Valdevit. Nanolattices: an emerging class of mechanical metamaterials. *Advanced Materials*, 29(40):1701850, 2017.

[3] A. Bensoussan, J.-L. Lions, and G. Papanicolaou. *Asymptotic Analysis for Periodic Structures*. AMS Chelsea Publishing, Providence, RI, 1978.

[4] A. Beskok and G. Em. Karniadakis. Simulation of heat and momentum transfer in complex microgeometries. *J. Thermal Physics & Heat Transfer*, 8(4):647–655, 1994.

[5] C.I. Christov. On a higher-gradient generalization of Fourier's law of heat conduction. In *AIP Conference Proceedings*, volume 946(1), pages 11–22. American Institute of Physics, 2007.
386

387

[6] B.D. Coleman and V.J. Mizel. Thermodynamics and departures from Fourier's law of heat conduction. *Arch. Rat. Mech. Anal.*, 13(1):245-261, 1963.

388

[7] W. Dahmen. On multivariate B-Splines. *SIAM J. Numer. Anal.*, 17(2):179-191, 1980.

389

[8] M. Gad el Hak. The fluid mechanics of microdevices: The Freeman scholar lecture. *J. Fluids Eng.*, 121(1):5-33, 1999.

390

[9] A. Fehér and R. Kovács. On the evaluation of non-Fourier effects in heat pulse experiments. *Int. J. Eng. Sci.*, 169:103577, 2021.

391

[10] A. Fehér, N. Lukács, L. Somlai, T. Fodor, M. Szücs, T. Fülöp, P. Ván, and R. Kovács. Size Effects and Beyond-Fourier Heat Conduction in Room-Temperature Experiments. *J. Non-Equil. Thermody.*, 46(4):403-411, 2021.

392

393

[11] S. Forest and M. Amestoy. Hypertemperature in thermoelastic solids. *C. R. Méc.*, 336(4):347-353, 2008.

394

[12] S. Forest, F. Pradel, and K. Sab. Asymptotic analysis of heterogeneous Cosserat media. *Int. J. Solids Structures*, 38(26-27):4585-4608, 2001.

395

396

[13] T. Frenzel, M. Kadic, and M. Wegener. Three-dimensional mechanical metamaterials with a twist. *Science*, 358(6366):1072-1074, 2017.

397

398

[14] L.R.T. Gardner and G.A. Gardner. A two dimensional bi-cubic b-spline finite element: used in a study of mhd-duct flow. *Computer Methods in Applied Mechanics and Engineering*, 124(4):365-375, 1995.

399

400

[15] P. Germain. La méthode des puissances virtuelles en mécanique des milieux continus, I: Théorie du second gradient. *J. Mécanique*, 12(2):235-274, 1973.

401

402

[16] P. Germain. The method of virtual power in continuum mechanics. Part 2: Microstructure. *SIAM J. Appl. Math.*, 25(3):556-575, 1973.

403

404

[17] G. Giantesio, A. Girelli, C. Lonati, A. Marzocchi, A. Musesti, and B. Straughan. Thermal convection in a higher velocity gradient and higher temperature gradient fluid. *J. Math. Fluid Mech.*, 27(3):47, 2025.

405

406

[18] M.E. Gurtin. *An introduction to Continuum Mechanics*, volume 158. Academic press, 1982.

407

[19] M.E. Gurtin. Generalized Ginzburg-Landau and Cahn-Hilliard equations based on a microforce balance. *Phys. D*, 92(3-4):178-192, 1996.

408

409

[20] K. Höllig. *Finite Element Methods with B-Splines*. Society for Industrial and Applied Mathematics, 2003.

410

[21] M. Kadic, T. Frenzel, and M. Wegener. When size matters. *Nature Physics*, 14(1):8-9, 2018.

411

[22] M. Lazar. On gradient field theories: gradient magnetostatics & gradient elasticity. *Philos. Mag.*, 94(25):2840-2874, 2014.

412

[23] M. Lazar, E. Agiasofitou, and T. Böhlke. Mathematical modeling of the elastic properties of cubic crystals at small scales based on the toupin-Mindlin anisotropic first strain gradient elasticity. *Contin. Mech. Thermodyn.*, 34(1):107-136, 2022.

413

414

[24] R. Lipton and B. Vernescu. Variational methods, size effects and extremal microgeometries for elastic composites with imperfect interface. *Math. Models Meth. Appl. Sci.*, 5:1139-1173, 08 1995.

415

416

[25] R. Lipton and B. Vernescu. Composites with imperfect interface. *Proc. R. Soc. Lond. A.*, 452:329-358, 1997.

417

[26] D. Matthes, E.-M. Rott, G. Savaré, and A. Schlichting. A structure preserving discretization for the Derrida-Lebowitz-Speer-Spohn equation based on diffusive transport. *Numer. Math.*, pages 1-49, 2025.

418

419

[27] G. Maugin. The method of virtual power in continuum mechanics: application to coupled fields. *Acta Mech.*, 35(1-2):1-70, 1980.

420

421

[28] C. C. Mei and B. Vernescu. *Homogenization Methods for Multiscale Mechanics*. World Scientific, 2010.

422

423

[29] R.D. Mindlin. Micro-structure in linear elasticity. *Arch. Rat. Mech. Anal.*, 16:51-78, 1964.

424

[30] R.D. Mindlin. On the equations of elastic materials with micro-structure. *Int. J. Solids Structures*, 1(1):73-78, 1965.

425

[31] R.D. Mindlin and N.N. Eshel. On first strain-gradient theories in linear elasticity. *Int. J. Solids Structures*, 4(1):109-124, 1968.

426

[32] D. Molavitabrizi, S. Khakalo, R. Bengtsson, and S.M. Mousavi. Second-order homogenization of 3-D lattice materials towards strain gradient media: numerical modelling and experimental verification. *Cont. Mechanics & Thermodyn.*, pages 1-20, 2023.

427

428

429

[33] G. Nika. A gradient system for a higher-gradient generalization of Fourier's law of heat conduction. *Mod. Phys. Lett. B*, 37(11):2350011, 2023.

430

431

[34] G. Nika. Derivation of effective models from heterogeneous Cosserat media via periodic unfolding. *Ric. Mat.*, 73(1):381-406, 2024.

432

433

[35] G. Nika. Scale-size dependent multi-continuum homogenization of complex bodies. *Quart. Appl. Math.*, 83:361-388, 2025.

434

435

[36] G. Nika and A. Muntean. Hypertemperature effects in heterogeneous media and thermal flux at small-length scales. *Netw. Heterog. Media*, 18(3):1207-1225, 2023.

436

[37] G. Nika and A. Muntean. Effective medium theory for second-gradient elasticity with chirality. *Asym. Analysis*, 139(1-2):111-137, 2024.

437

438

[38] W. Noll and B. Coleman. The thermodynamics of elastic materials with heat conduction and viscosity. In W. Noll, editor, *The Foundations of Mechanics and Thermodynamics: Selected Papers*, pages 145-156. Springer, 1974.

439

440

[39] R.W. Ogden. *Non-linear elastic deformations*. Ellis Harwood Ltd., Chichester, England, 1984.

441

442

[40] L. S. Pan, D. Xu, J. Lou, and Q. Yao. A generalized heat conduction model in rarefied gas. *Europhysics Letters (EPL)*, 73(6):846-850, 03 2006.

443

[41] P.M. Prenter. *Splines and variational methods*. Wiley, New York, 1975.

444

[42] E. Sanchez-Palencia. *Non-homogeneous media and vibration theory*. Springer-Verlag Berlin Heidelberg, 1980.

445

446

[43] B. Straughan. *Thermal convection with a Cattaneo flux law*, volume 53 of *Adv. Mech. Math*. Cham: Birkhäuser, 2024.

447

[44] R. Toupin. Elastic materials with couple-stresses. *Arch. Rat. Mech. Anal.*, 11(1):385-414, 1962.

448

[45] R. Toupin. Theory of elasticity with couple-stress. *Arch. Rat. Mech. Anal.*, 17:85-112, 1964.

449

[46] C. Truesdell and R. Toupin. The classical field theories. In *Principles of classical mechanics and field theory/Prinzipien der Klassischen Mechanik und Feldtheorie*, pages 226-858. Springer, 1960.

450