# Error Estimators for the Small-Biot Lumped Approximation
for the Conduction Dunking Problem

Kento Kaneko Massachusetts Institute of Technology, Cambridge, MA 02139, USA.
<EMAIL_ADDRESS><EMAIL_ADDRESS>Claude Le Bris École des Ponts and INRIA,
Champs-sur-Marne, 77455 Marne La Vallée, France.
<EMAIL_ADDRESS>Anthony T Patera Massachusetts Institute of Technology,
Cambridge, MA 02139, USA.
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We consider the dunking problem: a solid body at uniform temperature
$T_{\text{i}}$ is placed in a environment characterized by farfield
temperature $T_{\infty}$ and spatially uniform time-independent heat transfer
coefficient. We permit heterogeneous material composition: spatially dependent
density, specific heat, and thermal conductivity. Mathematically, the problem
is described by a heat equation with Robin boundary conditions. The crucial
parameter is the Biot number — a nondimensional heat transfer (Robin)
coefficient; we consider the limit of small Biot number.

We introduce first-order and second-order asymptotic approximations (in Biot
number) for several quantities of interest, notably the spatial domain average
temperature as a function of time; the first-order approximation is simply the
standard engineering ‘lumped’ model. We then provide asymptotic error
estimates for the first-order and second-order approximations for small Biot
number, and also, for the first-order approximation, alternative strict bounds
valid for all Biot number. Companion numerical solutions of the heat equation
confirm the effectiveness of the error estimates for small Biot number.

The second-order approximation and the first-order and second-order error
estimates depend on several functional outputs associated to an elliptic
partial differential equation; the latter is derived from Biot-sensitivity
analysis of the heat equation eigenproblem in the limit of small Biot number.
Most important is $\phi$, the only functional output required for the first-
order error estimates; $\phi$ admits a simple physical interpretation in terms
of conduction length scale. We investigate the domain and property dependence
of $\phi$: most notably, we characterize spatial domains for which the
standard lumped-model error criterion — Biot number (based on volume-to-area
length scale) small — is deficient.

Acknowledgments

We thank Prof. Masayuki Yano for his many contributions to the adaptive
refinement finite element method software suite used extensively in this work.

This work is supported by ONR Grant N00014-21-1-2382, Grant Monitor Dr Reza
Malek-Madani. The research of the second author is partially supported by ONR
and EOARD.

_Keywords:_ heat transfer, small Biot, lumped approximation, error estimation

##  1 Introduction

The dunking problem is ubiquitous in heat transfer engineering education and
professional practice: a body characterized by spatial domain
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\Omega}}$
and boundary
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\partial\Omega}}$,
at initial uniform temperature
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T_{\text{i}}}}$,
is immersed at time
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{t=0}}$ in
an environment (fluid and enclosure) at farfield temperature
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T_{\infty}}}$
[1, 2, 3]. The temperature distribution of the body,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}$,
evolves over the time interval
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{0<t\leq{t_{\text{final}}}}}$.
Note dimensional quantities shall be coded in blue, and nondimensional
quantities (to be introduced later) in black. In this introduction, for
simplicity of exposition, we largely consider homogeneous body composition,
and hence uniform thermophysical properties — in particular volumetric
specific heat,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\rho c}}$,
and scalar thermal conductivity,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{k}}$;
however, in subsequent development, we shall treat the important case of
material heterogeneity (and hence non-uniform thermophysical properties),
particularly in the form of several isotropic materials. The thermal
environment is characterized, most simply, by an average heat transfer
coefficient,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{h}}$.
Finally, we associate to the body domain two length scales: a chosen extrinsic
length scale,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\ell}}$;
an intrinsic length scale,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{L}\equiv|\Omega|/|\partial\Omega|}}$,
where
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{|\Omega|}}$
and
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{|\partial\Omega|}}$
denote the volume and surface area of the body, respectively. We note that the
use of the notation “$\equiv$” in this work signifies “is defined as”
(equivalent to the “:=” notation).

The solution
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T(\cdot,t)}}$
to the dunking problem satisfies, for
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{t\in(0,{t_{\text{final}}}]}}\,$,

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\rho
c\,\partial_{t}T=\boldsymbol{\nabla}\cdot(k\boldsymbol{\nabla}T)\text{ in
}\Omega\,,}}$ (1)

subject to boundary condition

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{k\,\partial_{n}T+h(T-T_{\infty})=0\text{
on }\partial\Omega}}\,,$ (2)

and initial condition

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T(\cdot,0)=T_{\text{i}}\text{
in }\Omega}}\,.$ (3)

Here, $\partial_{t}$ and $\partial_{n}$ refer to differentiation with respect
to time and domain boundary outward normal direction, respectively.

We nondimensionalize the spatial coordinate by
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\ell}}$
and time by
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\rho
c\ell^{2}/k}}$. The nondimensional temperature is defined in terms of the
dimensional temperature
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}$ as
$u\equiv{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{(T-T_{\infty})/({T_{\text{i}}}-T_{\infty})}}$.
The mathematical formulation of the nondimensional dunking problem then
corresponds to a heat equation characterized by the spatial domain $\Omega$
and a single parameter, the Robin coefficient. The Robin coefficient is given
by an extrinsic Biot number,
$B\equiv{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{h\ell/k}}$;
we also define an intrinsic Biot number,
$\text{Bi}_{\text{dunk}}\equiv{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{h\mathscr{L}/k}}=B\mathscr{L}$,
where
$\mathscr{L}={\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{L}/\ell}}$
(note $\text{Bi}_{\text{dunk}}$ is typically indicated in textbooks simply by
Bi [1, 2, 3]). As we shall see, the general case of heterogeneous properties
is considerably more complicated: the nondimensional volumetric specific heat
and thermal conductivity will now be functions of the spatial coordinate.

In this work, we study the limit of small Biot number. The small-Biot limit
arises quite often in practice in particular for “everyday” artifacts at
modest temperatures subject to natural convection and radiation. In contrast,
larger systems at very elevated temperatures subject to brisk forced
convection or change-of-phase typically yield larger Biot numbers. An example
of a small-Biot application is the process of annealing by natural convection
— for instance, to relieve residual stresses.

We consider two quantities of interest: the (spatial) domain average
temperature as a function of time,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T_{\text{avg}}(t)}}$
— also directly related to the heat loss (or gain) of the body from (or to)
the environment; the (spatial) domain-boundary average temperature as a
function of time,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T_{\partial{\text{avg}}}(t)}}$.
The nondimensional versions of these quantities of interest are given
respectively by
$u_{\text{avg}}(t)\equiv{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{(T_{\text{avg}}(}}t{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\rho
c\ell^{2}/k)-T_{\infty})/(T_{\text{i}}-T_{\infty})}}$ and
$u_{\partial{\text{avg}}}(t)\equiv{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{(T_{\partial{\text{avg}}}(}}t{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\rho
c\ell^{2}/k)-T_{\infty})/(T_{\text{i}}-T_{\infty})}}$; it shall also prove
convenient to define
$u_{\Delta}(t)\equiv[({u_{\text{avg}}}-{u_{\partial{\text{avg}}}})/{u_{\text{avg}}}](t)$,
which measure the relative temperature variation within the body.

The first-order111Note in general “order” refers to the convergence rate in
$B$ (or, equivalently, $\text{Bi}_{\text{dunk}}$). ‘classic’ lumped
approximation to $u_{\text{avg}}$ for the small-$\text{Bi}_{\text{dunk}}$
dunking problem is a simple exponential in time,
$\tilde{u}_{\text{avg}}^{1}(t)\equiv\exp(-B\gamma t)$ [1, 2, 3], where
$\gamma\equiv 1/\mathscr{L}$. However, and despite the simplicity and utility
of this small-$\text{Bi}_{\text{dunk}}$ lumped approximation, there is
relatively little analysis of the associated approximation error. For example,
most textbooks [1, 2, 3] indicate only that $\text{Bi}_{\text{dunk}}$ must be
small; in some cases, a threshold might be provided, for example
$\text{Bi}_{\text{dunk}}\leq 0.1$ — but typically qualified by “usually”. A
notable exception is the important work of Gockenbach and Schmidtke [4], in
which a rigorous asymptotic error estimate is provided for the particular case
of a sphere with homogeneous thermophysical properties and uniform heat
transfer coefficient:
$|\tilde{u}_{\text{avg}}^{1}(t)-u_{\text{avg}}(t)|\leq(3/5)\,\text{Bi}_{\text{dunk}}/e$
for all time $t$. The analysis of [4] is based on an eigenvalue expansion in
Biot.

We enumerate here the specific contributions of this work:

  1. 1.

We consider small-$\text{Bi}_{\text{dunk}}$ approximation and error estimation
based on asymptotic expansion of the first eigenfunction (for the eigenproblem
associated with the heat equation),
$\psi_{1}(\cdot;B)=\psi_{1}(\cdot;0)+B\frac{d\psi_{1}}{dB}(\cdot;0)+\pazocal{O}(B^{2})$,
where $\frac{d\psi_{1}}{dB}(\cdot;0)$ is found following standard sensitivity-
derivative approaches [5, 6] as the solution to an elliptic partial
differential equation (PDE). The key quantity in our analysis is
$\phi\equiv\int_{\Omega}\kappa\nabla\frac{d\psi_{1}}{dB}(\cdot;0)\cdot\nabla\frac{d\psi_{1}}{dB}(\cdot;0)$,
where $\kappa$ is the nondimensionalized thermal conductivity based on the
infimum of
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{k}}$ over
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\Omega}}$.
For canonical domains and selected triangles, we provide closed-form solutions
for $\phi$; we also provide an efficient computational framework for
evaluation of $\phi$ for general domains.

  2. 2.

We introduce a new second-order lumped approximation to the solution of the
thermal dunking problem, $\tilde{u}_{\text{avg}}^{2\text{P}}$, which is
dependent on $\phi$. The new approximation has improved accuracy over the
‘classic’ approximation, $\tilde{u}_{\text{avg}}^{1}$, for general geometry
and general thermophysical property distribution. We also introduce an
approximation for the spatial temperature variation, $\tilde{u}_{\Delta}^{2}$.

  3. 3.

We present an error analysis for all three approximations:
$\tilde{u}_{\text{avg}}^{1}$, $\tilde{u}_{\text{avg}}^{2\text{P}}$, and
$\tilde{u}_{\Delta}^{2}$. The error analysis can be performed at marginal
computational cost over the cost of computing $\phi$. For the first-order
approximation, we can provide two additional results: we demonstrate that
$\tilde{u}_{\text{avg}}^{1}$ is in fact a lower bound of $u_{\text{avg}}$ for
all admissible $t$ and $B$; we provide a strict error bound for
$\tilde{u}_{\text{avg}}^{1}$ which is valid for all $B$.

  4. 4.

We provide an explicit engineering interpretation of $\phi$ in terms of an
improved engineering conduction length scale associated with a body thermal
resistance. The latter generalizes in a rigorous way the usual definition of
thermal resistance — between two surfaces — to the notion of thermal
resistance between body effective center and body boundary. Within this
interpretation framework, we motivate $\tilde{u}_{\text{avg}}^{1}$,
$\tilde{u}_{\text{avg}}^{2\text{P}}$, and $\tilde{u}_{\Delta}^{2}$.

  5. 5.

We provide an extensive discussion of the dependence of $\phi$ on spatial
domain — in fact, only body shape matters — and body composition. In some
cases we can provide useful approximations and bounds for $\phi$ which require
relatively little information about (say) the spatial distribution of the
thermophysical properties,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\rho c}}$
and ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{k}}$.
We define a distance between spatial domains which predicts well the stability
of $\phi$ with respect to geometric perturbation. More generally, we identify
two geometry features in terms of which we can characterize spatial domains
for which $\phi$ is large — and the textbook lumped error criterion
potentially misleading.

  6. 6.

We provide a numerical study of $\phi$ with emphasis on geometry and
thermophysical properties. The numerical results are based on an adaptive
finite element framework [7] with associated a posteriori error estimators;
the latter ensure that the finite element error does not compromise any of our
conclusions. We confirm (i) the effectiveness of our error estimators, and
(ii) the relevance of our distance and geometry features in understanding the
dependence of $\phi$ on spatial domain.

We provide more extensive references in context.

There is one aspect in which our analysis is incomplete: the treatment of the
heat transfer coefficient,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{h}}$. We
consider here almost exclusively the case in which the heat transfer
coefficient is independent of time and furthermore uniform in space. For the
case of natural convection and radiation, the heat-transfer coefficient is
implicitly linearized about the initial temperature of the body. The more
realistic case of temporally-dependent and spatially-dependent heat transfer
coefficient is best addressed in the context of a more general discussion
which takes as point of departure the “truth” conjugate heat transfer
formulation; this approach will be reported in a subsequent publication.
However, we provide even here, in Appendix D, a simple first result.

We now provide a roadmap for the manuscript.

  * In Section 2 we present the $\text{Bi}_{\text{dunk}}$ initial-value-problem formulation as well as the associated $\text{Bi}_{\text{dunk}}$ eigenproblem; we also develop the connection between the former and the latter through a separation of variables expansion.

  * In Section 3 we then develop the small-Biot asymptotic analysis of the eigenproblem: as $B\rightarrow 0$, (i) $\lambda_{1}$ (the first eigenvalue) $\sim B\gamma-\phi B^{2}$, and (ii) the spatial structure is well-approximated by $\psi_{1}$ (the first eigenmode). We also provide in Section 3 the sensitivity formulation for $\phi$. In Section 5 we then analyze the sensitivity formulation for $\phi$ as regards general properties and particular illustrative solutions. Note many of the results of Section 3.1 are particular instantiations of more general eigenproblem perturbation techniques [5, 6].

  * In Section 4 we develop, through the resistance formulation, the lumped approximations $\tilde{u}_{\text{avg}}^{1}(t)$, $\tilde{u}_{\text{avg}}^{2}(t)$, and $\tilde{u}_{\Delta}^{2}$. In Section 6 we then analyze the error in these lumped approximations and provide, for all approximations, asymptotic error estimates (as $B\rightarrow 0$) as well as strict error bounds, valid for all $B$, for $\tilde{u}_{\text{avg}}^{1}$.

Note we intercalate the discussion of $\phi$ and the lumped approximation for
pedagogical purposes.

The main part of the manuscript, intended for a more general heat transfer
audience, emphasizes the presentation and interpretation of the error
estimates; we include the statement of the theoretical results but without
proof. The appendices to the manuscript include all necessary mathematical
details as well as sketches of the proofs of the main propositions. We note
that many of the mathematical proofs are indeed elementary and many of the
results are well known, however we provide them here for completeness and for
the convenience of the reader.

##  2 Governing Equations

###  2.1 Initial Value Problem

A passive solid body at uniform temperature
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T_{\text{i}}}}$
is abruptly placed — dunked — at time
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{t=0}}$ in
an environment, fluid and enclosure, at initial and farfield temperature
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T_{\infty}}}$.
Note that blue variables are dimensional, and black variables are
nondimensional. Most of the analysis presented in this work can be extended to
more general contexts, for example uniform heat generation in the body,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{q^{\textsc{v}}(t)}}$,
or time-dependent farfield temperature,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T_{\infty}(t)}}$.
Note that, in general, units are SI, and temperature is in degrees Celsius (or
Kelvin).

We first introduce the spatial domain of the body we consider,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\Omega}}\subset\mathbb{R}^{d}$,
and associated boundary
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\partial\Omega}}$;
we consider $d\in\\{1,2,3\\}$. We shall require
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\Omega}}$
Lipschitz. A point in
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\Omega}}$
shall be denoted
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{x\equiv(x_{1},\ldots,x_{d})}}$.
We choose for our length scale some characteristic dimension of
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\Omega}}$,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\ell}}$,
for example the diameter or the InRadius. We also introduce an intrinsic
length scale
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{L}\equiv|\Omega|/|\partial\Omega|}}$,
and an intrinsic associated inverse length scale,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\gamma\equiv|\partial\Omega|/|\Omega|}}$;
here $|\,\cdot|$ denotes measure, and hence, in dimension $d=3$,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\gamma}}$
is the dimensional surface area to volume ratio. We shall denote time by
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{t}}$; we
shall restrict attention to the temporal interval of interest
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{t\in[0,{t_{\text{final}}}]}}$.

We next introduce two nondimensional quantities related to thermophysical
properties:

$\displaystyle\sigma=\dfrac{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\rho
c}}}{\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\Omega}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{(\rho
c)}}}\,,$ (4)
$\displaystyle\kappa={\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\dfrac{k}{k_{\inf}}}}\,,$
(5)

where ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\rho
c}}$ and
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{k}}$ are,
respectively, the volumetric specific heat and thermal conductivity, both
strictly positive and scalar-valued (generalization to anisotropic heat-
transfer with symmetric positive definite second-order $k$ tensor is not a
subject of this work). Note
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\rho}}$ is
the density and
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{c}}$ is
the mass specific heat; the product
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\rho c}}$
is the volumetric specific heat (which is bounded away from zero). Here
$\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int$ refers to the average of the integrand over the
domain of integration, and
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{k_{\inf}}}$
(more properly
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{k_{\operatorname*{ess\,inf}}}}$)
is the infimum of
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{k}}$ over
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\Omega}}$.
It follows from our definitions that $\sigma>0$ and
$\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\Omega}\sigma=1$, (and $\kappa>0$) and
$\inf_{\Omega}\kappa=1$. Note also that $\sigma=1$ corresponds to spatially
uniform volumetric specific heat, and $\kappa=1$ corresponds to spatially
uniform thermal conductivity.

We may then define the extrinsic Biot number as

$\displaystyle
B\equiv{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\frac{h\ell}{k_{\inf}}}}\,,$
(6)

and the intrinsic Biot number as

$\displaystyle\mathrm{Bi}_{\mathrm{dunk}}={\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\dfrac{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{h\mathscr{L}}}}{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{k_{\inf}}}}}}\;;$
(7)

and hence
$\mathrm{Bi}_{\mathrm{dunk}}=B({\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{L}/\ell}})=B/(\ell{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\gamma}})$.
Note that

$\displaystyle\gamma\equiv{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\ell\gamma}}$
(8)

shall appear prominently in the exposition. As indicated in the introduction,
we consider in the present work, the case in which the heat transfer
coefficient
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{h}}$ is
taken as uniform in space and constant in time. We further assume that the
heat transfer coefficient is non-negative, and hence also $B$ is non-negative.

We now define the nondimensional spatial coordinate as
$x\equiv{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{x/\ell}}$
and the nondimensional temporal coordinate as
$t\equiv{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{t/t_{\text{diff}}}}$,
where

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{t_{\text{diff}}\equiv}}\dfrac{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\ell^{2}}}\,\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\Omega}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{(\rho
c)}}}{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{k_{\inf}}}}\,.$
(9)

We further introduce the nondimensional temperature

$\displaystyle
u\equiv{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\dfrac{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}-{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T_{\infty}}}}{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T_{\text{i}}}}-{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T_{\infty}}}}}}\,.$
(10)

Our nondimensionalization reduces to the standard textbook form in the case in
which the thermophysical properties are uniform; notably, $t$ is the Fourier
number,
$\text{Fo}_{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\ell}}}$.

The nondimensional temperature, $u(x,t;B;\sigma,\kappa)$, then satisfies, for
any $B\in\mathbb{R}_{0+}\equiv\\{x\in\mathbb{R}\,|\,x\geq 0\\}$,

$\displaystyle\sigma\partial_{t}u=\boldsymbol{\nabla}\cdot(\kappa\boldsymbol{\nabla}u)\text{
in }\Omega\,,\;0<t\leq{t_{\text{final}}}\,,$ (11)

subject to boundary condition

$\displaystyle\kappa\partial_{n}u+Bu=0\text{ on
}\partial\Omega\,,\;0<t\leq{t_{\text{final}}}\,,$ (12)

and initial condition

$\displaystyle u=1\text{ in }\Omega\text{ at }t=0\;.$ (13)

Here $\partial_{t}$ and $\partial_{n}$ refer to the partial derivatives with
respect to time and the outward normal direction, respectively. We note that,
strictly speaking, the strong formulation (11) – (13) requires $\sigma$ and
$\kappa$ sufficiently smooth. However, the corresponding weak formulation,
provided in Appendix A.2, admits discontinuous data, and hence our analysis —
based on the weak form — shall directly apply to the practically important
case of material heterogeneity.

We shall in general suppress the $\sigma$ and $\kappa$ parameters unless we
are specifically focused on the dependence on these variables. We shall write
$u(\cdot,t;B)$ or simply $u(t;B)$ to indicate the spatial nondimensional
temperature function for any time $t$ and Biot number $B$ (and, implicitly,
prescribed $\sigma$ and $\kappa$); we shall write $u(x,t;B)$ to refer to the
evaluation of the spatial nondimensional temperature function at spatial
coordinate $x\in\Omega$.

###  2.2 Eigenproblem

We can associate to our initial value problem, in standard fashion, the
corresponding eigenproblem: for given $B\in\mathbb{R}_{0+}$, find
$(\psi_{j},\lambda_{j})_{j=1,2,\ldots}$ solution to

$\displaystyle-\boldsymbol{\nabla}\cdot(\kappa\boldsymbol{\nabla}\psi_{j})=\lambda_{j}\sigma\psi_{j}\text{
in }\Omega\,,$ (14)

subject to boundary condition

$\displaystyle\kappa\partial_{n}\psi_{j}+B\psi_{j}=0\text{ on
}\partial\Omega\,,$ (15)

and normalization

$\displaystyle\int_{\Omega}\sigma\psi^{2}_{j}=1\;.$ (16)

The corresponding weak form is provided in Appendix A.3. We shall write
$\psi_{j}(\cdot;B)$ or simply $\psi_{j}(B)$ to indicate the spatial
eigenfunction for any given Biot number $B$ (and implicitly, prescribed
$\sigma$ and $\kappa$); we shall write $\psi_{j}(x;B)$ to refer to the
evaluation of the spatial eigenfunction function at spatial coordinate
$x\in\Omega$.

It is a standard result that the eigenvalues are non-negative for $B\geq 0$.
We shall enumerate the eigenvalues in increasing magnitude, hence
$0\leq\lambda_{1}(B)<\lambda_{2}(B)\leq\lambda_{3}(B)\leq\ldots$. It is also
standard to note that the eigenfunctions are orthonormal:

$\displaystyle\int_{\Omega}\sigma\psi_{i}(B)\,\psi_{j}(B)=\delta_{ij}\;,\;i,j=1,2,\dots\;;$
(17)

here $\delta_{ij}$ is the Kronecker-delta symbol.

Considering $\lambda_{1}$, it can be shown that $\lambda_{1}$ is of
multiplicity unity, and that (i) $\psi_{1}$ has constant sign, (ii) $\psi_{1}$
is (up to a change of sign) strictly positive, and (iii) $\psi_{1}$ is (up to
a change of sign) unique. As regards (ii) and (iii), see for example [8]. We
shall choose $\psi_{1}(B)$ positive, which, with (16), uniquely specifies
$\psi_{1}(B)$. We shall henceforth adopt the convention that $\psi(B)$
(respectively, $\lambda(B)$) without subscript shall refer to the first
eigenfunction, $\psi_{1}(B)$ (respectively, the first eigenvalue,
$\lambda_{1}(B)$), and that furthermore $\psi^{0}$ (respectively,
$\lambda^{0}$) shall refer to $\psi(B=0)$ (respectively, $\lambda(B=0))$.

###  2.3 Separation of Variables Representation

Much of our analysis shall rely on the classical separation of variables
representation of our solution. We begin with a lemma which shall serve
subsequently in several contexts:

######  Lemma 2.1 (Representation of Unity).

We may write, for any $B\in\mathbb{R}_{0+}$,

$\displaystyle\lim_{N\rightarrow\infty}\,\left|\left|\;1-\sum_{k=1}^{N}|\Omega|M(\psi_{k}(B))\psi_{k}(\cdot;B)\;\right|\right|_{L^{2}(\Omega)}\\!\rightarrow
0\;,$ (18)

or, in abbreviated form (to be interpreted as (18))

$\displaystyle
1=\sum_{k=1}^{\infty}|\Omega|M(\psi_{k}(B))\psi_{k}(\cdot;B)\;;$ (19)

here $\|w\|_{L^{2}(\Omega)}\equiv(\int_{\Omega}w^{2})^{1/2}$. Furthermore,

$\displaystyle\lim_{N\rightarrow\infty}\,|\Omega|\sum_{k=1}^{N}(M(\psi_{k}(B)))^{2}=1\;,$
(20)

or, in abbreviated form (to be interpreted as (20))

$\displaystyle 1=|\Omega|\sum_{k=1}^{\infty}(M(\psi_{k}(B)))^{2}\;.$ (21)

Here $|\Omega|$ is the measure of $\Omega$ and $M:X\rightarrow\mathbb{R}$ is
given by

$\displaystyle M(v;\sigma)\equiv\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\Omega}\sigma v\;.$ (22)

We note that for $\sigma=1$ (uniform volumetric specific heat) $M$ evaluates
to the domain mean of $v$. Recall also that $\psi(B)$ refers to the function
$\psi(\cdot;B)\in X$, whereas $\psi(x;B)$ refers to the evaluation of
$\psi(\cdot;B)$ at spatial coordinate $x\in\Omega$.

Hence $X\equiv H^{1}(\Omega)$, where $H^{1}(\Omega)$ is imbued with norm
$||\cdot||_{H^{1}}\equiv||\cdot||_{X}=\left(\int_{\Omega}\cdot\right)^{1/2}$

The sketch of the proof of Lemma 2.1 is provided in Appendix C.1.

We can then write our solution as a standard separation of variables sum:

######  Proposition 2.2 (Separation of Variables: Field).

For any $B\in\mathbb{R}_{0+}$,

$\displaystyle
u^{N}(\cdot,t;B)\equiv\sum_{k=1}^{N}|\Omega|M(\psi_{k}(B))\psi_{k}(\cdot;B)\exp(-\lambda_{k}(B)t)\;.$
(23)

Then, for all $t\in(0,{t_{\text{final}}}]$,

$\displaystyle\lim_{N\rightarrow\infty}\|\;u^{N}(\cdot,t;B)-u(\cdot,t;B)\;\|_{X}\rightarrow
0\;,$ (24)

or in abbreviated form (to be interpreted as (24))

$\displaystyle
u(\cdot,t;B)\equiv\sum_{k=1}^{\infty}|\Omega|M(\psi_{k}(B))\psi_{k}(\cdot;B)\exp(-\lambda_{k}(B)t)\;.$
(25)

Recall that the $X$ norm is the $H^{1}(\Omega)$ norm.

The sketch of the proof of Proposition 2.2 is provided in Appendix C.2. We
note here only that $u(\cdot,t;B)$ of (25) for $t=0$ reduces to
$\sum_{k=1}^{\infty}|\Omega|M(\psi_{k}(B))\psi_{k}(\cdot;B)=1$ from Lemma 2.1.

###  2.4 Quantities of Interest

Quantity of Interest (QoI) shall refer to a functional of the temperature
$u(\cdot,t;B)$. We shall consider in the present work three QoI: the domain
mean, denoted $u_{\text{avg}}$; the boundary mean, denoted
$u_{\partial{\text{avg}}}$; and the relative domain mean-boundary mean
difference, $u_{\Delta}$. In our analysis and results, we focus on
$u_{\text{avg}}$ and $u_{\Delta}$; $u_{\partial{\text{avg}}}$ is a necessary
intermediate definition. The domain mean QoI is defined as

$\displaystyle u_{\text{avg}}(t;B)\equiv
M(u(\cdot,t;B))=\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\Omega}\sigma u(\cdot,t;B)\;;$ (26)

the boundary mean QoI is defined as

$\displaystyle u_{\partial{\text{avg}}}(t;B)\equiv
H(u(\cdot,t;B))=\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\partial\Omega}u(\cdot,t;B)\;;$ (27)

the relative domain mean-boundary mean difference QoI is defined as

$\displaystyle
u_{\Delta}(t;B)\equiv\dfrac{u_{\text{avg}}(t;B)-u_{\partial{\text{avg}}}(t;B)}{u_{\text{avg}}(t;B)}\;.$
(28)

Here, $M$ is defined in (22), and $H:X\rightarrow\mathbb{R}$ is given by

$\displaystyle H(v)\equiv\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\partial\Omega}v\;;$ (29)

note $H$ evaluates to the boundary mean of $v$. These QoIs are defined for
almost all times $t$ in our interval $[0,{t_{\text{final}}}]$. Note it follows
from the maximum principle [9] that $u_{\text{avg}}$ and
$u_{\partial{\text{avg}}}$ take on values in the range $[0,1]$.

We can now express our QoI in terms of our separation of variables
representation:

######  Proposition 2.3 (Separation of Variables: QoI).

We may write

$\displaystyle
u_{\text{avg}}(t;B)=\sum_{k=1}^{\infty}|\Omega|(M(\psi_{k}(B)))^{2}\exp(-\lambda_{k}(B)t)\;$
(30)

and

$\displaystyle
u_{\partial{\text{avg}}}(t;B)=\sum_{k=1}^{\infty}|\Omega|\,M(\psi_{k})H(\psi_{k})\exp(-\lambda_{k}(B)t)$
(31)

or

$\displaystyle
u_{\partial{\text{avg}}}(t;B)=\sum_{k=1}^{\infty}|\Omega|(M(\psi_{k}(B)))^{2}\,\dfrac{\lambda_{k}(B)}{B\gamma}\exp(-\lambda_{k}(B)t)\;,$
(32)

where $M$ and $H$ are defined in (22) and (29), respectively; recall also that
$\gamma\equiv|\partial\Omega|/|\Omega|$ is the (nondimensional) surface area
to volume ratio as given in (8).

We provide a sketch of the proof in Appendix C.3.

We note that the QoI, in particular the domain-mean QoI, can serve in either
forward or inverse mode. In forward mode, given $B$, we wish to evaluate
$u_{\text{avg}}(t;B),0\leq t\leq{t_{\text{final}}}$. In inverse mode, given a
target domain-mean QoI value, $u^{*}_{\text{avg}}$, we wish to evaluate
$\tau_{\text{avg}}(B)$ such that

$\displaystyle
u_{\text{avg}}(\tau_{\text{avg}}(u_{\text{avg}}^{*};B);B)=u_{\text{avg}}^{*}\;.$
(33)

For example, if $u_{\text{avg}}^{*}=\exp(-1)$, then we might interpret
$\tau_{\text{avg}}$ as the time constant of the system. In some cases
$\tau_{\text{avg}}$ can be a more meaningful and sensitive output than
$u_{\text{avg}}$. The error estimators for $\tau_{\text{avg}}$ follow directly
from the error estimators for $u_{\text{avg}}$; we shall not further consider
the inverse mode in this work.

##  3 Eigenproblem for Small $B$

###  3.1 Asymptotic Analysis

We first introduce some notation. As before, $\psi$ and $\lambda$ shall refer
respectively to $\psi_{1}$ and $\lambda_{1}$, and $\psi^{0}$ and $\lambda^{0}$
shall refer respectively to $\psi_{1}$ and $\lambda_{1}$ evaluated at $B=0$.
We now further define $\psi^{\prime}$ as the first derivative of $\psi_{1}$
with respect to $B$, $d\psi/dB(B)$; similarly, we define $\lambda^{\prime}$ as
the first derivative of $\lambda_{1}$ with respect to $B$, $d\lambda/dB(B)$.
We define $\psi^{\prime\,0}$ as the first derivative of $\psi_{1}$ with
respect to $B$ evaluated at $B=0$, $\psi^{\prime}(B=0)$; similarly, we define
$\lambda^{\prime\,0}$ as the first derivative of $\lambda_{1}$ with respect to
$B$ evaluated at $B=0$, $\lambda^{\prime}(B=0)$. Finally, second derivatives
with respect to $B$ are denoted by ′′ and also inherit 0 for evaluation at
$B=0$. The existence of the derivatives of $\psi_{1}$ and $\lambda_{1}$ with
respect to $B$, in fact to all orders, follows from the uniqueness of the
first eigenfunction and continuity considerations.

As we shall perturb about $B=0$, we first summarize the case $B=0$:

######  Lemma 3.1 ((Eigenfunction, Eigenvalue) for $B=0$).

For $B=0$ we obtain $\psi^{0}(\equiv\psi_{1}(B=0))=|\Omega|^{-1/2}$ and
$\lambda^{0}(\equiv\lambda_{1}(B=0))=0$.

Note that for $B=0$ the Robin problem reduces to the Neumann problem. The
proof of Lemma 3.1 is provided in Appendix C.4.

The small-$B$ analysis will depend on the predominance of the first mode. In
preparation, we provide the following lemma:

######  Lemma 3.2 (Amplitude of First Eigenfunction).

For $B\rightarrow 0$, the square of the first eigenfunction
($\sigma$–weighted) domain mean admits expansion

$\displaystyle|\Omega|(M(\psi(B)))^{2}=1-B^{2}\Upsilon+\pazocal{O}(B^{3})\;,$
(34)

where

$\displaystyle\Upsilon\equiv\int_{\Omega}\sigma\,\psi^{\prime\,0}\,\psi^{\prime\,0}\;$
(35)

is non-negative and independent of $B$.

We may thus conclude that, indeed, as $B\rightarrow 0$, the spatial structure
is well-approximated by the first mode. We provide a sketch of the proof in
Appendix C.5.

We now provide the main result of this section:

######  Proposition 3.3 (First Eigenvalue Asymptotic Expansion).

The first eigenvalue admits expansion

$\displaystyle\lambda(B;\kappa,\sigma)=B\gamma-\phi
B^{2}+\pazocal{O}(B^{3})\text{ as }B\rightarrow 0\;,$ (36)

where $\phi\in\mathbb{R}_{0+}$ is given by

$\displaystyle\phi(\sigma,\kappa)\equiv\int_{\Omega}\kappa\boldsymbol{\nabla}\psi^{\prime\,0}(\cdot;\sigma,\kappa)\cdot\boldsymbol{\nabla}\psi^{\prime\,0}(\cdot;\sigma,\kappa)\;;$
(37)

we recall the definition $\psi^{\prime\,0}\equiv d\psi_{1}/dB\text{ evaluated
at }B=0$.

We provide in Appendix C.6 a sketch of the proof of Proposition 3.3.

We can then readily develop a first-order approximation to
$\lambda(B)(\equiv\lambda(B;\kappa,\sigma))$ as $B\rightarrow 0$:

######  Proposition 3.4 (First Eigenvalue: First-Order and Second-Order
Approximations in $B$).

We introduce the approximations

$\displaystyle\tilde{\lambda}^{1}(B)$ $\displaystyle\equiv B\gamma\;,$ (38)
$\displaystyle\tilde{\lambda}^{2}(B)$ $\displaystyle\equiv B\gamma-\phi
B^{2}\;,$ (39) $\displaystyle\tilde{\lambda}^{2\text{P}}(B)$
$\displaystyle\equiv\frac{B\gamma}{1+\phi B/\gamma}\;.$ (40)

Here, ‘P’ in $\tilde{\lambda}^{2\text{P}}(B)$ refers to Padé as
$\tilde{\lambda}^{2\text{P}}(B)$ is the $[1/1]$ Padé approximant [10] to
$\lambda(B)$. These approximations satisfy

$\displaystyle\left|\dfrac{\lambda(B)-\tilde{\lambda}^{1}(B)}{\lambda(B)}\right|=\phi\underbrace{(B/\gamma)}_{\text{Bi}_{\text{dunk}}}+\pazocal{O}(B^{2})\text{
as }B\rightarrow 0\;,$ (41)
$\displaystyle\left|\dfrac{\lambda(B)-\tilde{\lambda}^{2}(B)}{\lambda(B)}\right|=\pazocal{O}(B^{2})\text{
as }B\rightarrow 0\;,$ (42)
$\displaystyle\left|\dfrac{\lambda(B)-\tilde{\lambda}^{2\text{P}}(B)}{\lambda(B)}\right|=\pazocal{O}(B^{2})\text{
as }B\rightarrow 0\;,$ (43)

where we recall that $\text{Bi}_{\text{dunk}}\equiv B/\gamma$ is the standard
engineering definition of the Biot number in the dunking context [1, 2, 3].
The superscript numerals in $\tilde{\lambda}^{1}$, $\tilde{\lambda}^{2}$, and
$\tilde{\lambda}^{2\text{P}}$ refer to the leading-order of the relative
approximation error.

The proof of Proposition 3.4 follows directly from Proposition 3.3: evaluating
the Taylor series around $B=0$ of the left-hand sides of (41) – (43).

###  3.2 Sensitivity Formulation for $\phi$

We observe that $\phi$ plays a key role in the error estimate for
$\tilde{\lambda}^{1}$ — and, in the next section, the error estimate for the
lumped approximation to the initial-value problem — and also in the definition
of $\tilde{\lambda}^{2\text{P}}$. In turn, $\phi$ is determined by the
function $\psi^{\prime\,0}$, which we shall denote the sensitivity derivative.
In fact, the definition of $\psi^{\prime\,0}$ as the derivative of the (first)
eigenfunction at $B=0$ is not directly amenable to either easy analysis or
efficient/stable computation. We thus present an alternative formulation:

######  Proposition 3.5 (Elliptic Equation for Sensitivity Derivative).

The sensitivity derivative $\psi^{\prime\,0}$ is the solution to

$\displaystyle-\boldsymbol{\nabla}\cdot(\kappa\boldsymbol{\nabla}\psi^{\prime\,0})=|\Omega|^{-1/2}\gamma\sigma\text{
in }\Omega\;$ (44)

subject to boundary condition

$\displaystyle\kappa\partial_{n}\psi^{\prime\,0}=-|\Omega|^{-1/2}\text{ on
}\partial\Omega\;$ (45)

and zero-mean condition

$\displaystyle\int_{\Omega}\sigma\psi^{\prime\,0}=0\;.$ (46)

We note that (44) – (46) is solvable — $\sigma$ is of mean unity over $\Omega$
— and admits a unique solution. We recall that $\partial_{n}$ denotes the
derivative with respect to the outward normal.

We shall refer to (44) – (46) as the (strong form of the) “Sensitivity
Equation”, or more elaborately “Equation for Sensitivity of the Eigenfunction
to Biot”. We provide in Appendix A.4 the weak form of the Sensitivity
Equation, and in Appendix C.7 a sketch of the proof of Proposition 3.5.

##  4 Lumped Formulation

In this section we apply engineering thermal resistance concepts to (i)
understand the physical significance of $\phi$, and (ii) develop our lumped
approximations. Note in this section there is no pretense of rigor: error
estimates for the lumped approximations — which serve as justifications for
our choices — are provided in Section 6 . In Section 4.2, we shall provide the
approximate quantities introduced in the derivations with $\hat{\cdot}$ to
emphasize the formal nature of the discussion. The main result of this section
is the engineering motivation for the three quantities of interest: (i) the
(first-order) ‘classical’ lumped approximation to $u_{\text{avg}}$, (ii) a
second-order lumped approximation to $u_{\text{avg}}$, and (iii) an
approximation to $u_{\Delta}$ which is null under the (first-order)
‘classical’ theory. These approximations are summarized in Definition 4.2.

In this section, we first introduce the body average resistance
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{R}_{\text{avg}}}}$:
the effective resistance between the body average temperature and the surface
temperature. In this resistance analogy, we derive an effective conduction
length scale $\mathscr{L}_{\text{cond}}$ based on the standard intrinsic
length scale $\mathscr{L}$ but now corrected for second-order effects. We then
move to a thermal circuit analogy to derive the effective (corrected) Biot
number $\text{Bi}^{\prime}_{\text{dunk}}$. Note both corrections are based on
$\phi$. Finally, these corrected quantities are then incorporated in a second-
order approximation of the solution, which improves upon the (first-order)
‘classic’ approximation.

###  4.1 Body Average Resistance

We first introduce a conduction problem in which (non-uniform) heat generation
balances (uniform) heat flux at the boundary:

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{-\boldsymbol{\nabla}\cdot(k\boldsymbol{\nabla}T^{*})}}\
$
$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{=}}\
\sigma{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{q\gamma}}\text{
in
}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\Omega}},$
(47)
$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{k\partial_{n}T^{*}}}\
$
$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{=-q}}\text{
on
}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\partial\Omega}},$
(48) $\displaystyle\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\Omega}\sigma{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{(T^{*}-T_{\text{ref}})}}\
$
$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{=0}}\,;$
(49)

here
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T_{\text{ref}}}}$
is a reference temperature. We can demonstrate that equations (47) – (49) are
solvable: we integrate (47) over
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\Omega}}$,
apply (48), and recall the definition of
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\gamma}}$
(8) of Section 2. We now define the associated body average thermal resistance
as

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{R}_{\text{avg}}}}$
$\displaystyle\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\equiv}}\
\dfrac{\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\Omega}\sigma{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T^{*}-}}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\partial\Omega}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T^{*}}}}{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{q|\partial\Omega|}}}$
(50) $\displaystyle\
{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{=}}\
\dfrac{\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\Omega}\sigma{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{(T^{*}-T_{\text{ref}})}}-\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\partial\Omega}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{(T^{*}-T_{\text{ref}})}}}{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{q|\partial\Omega|}}},$
(51)

in which the second equality (51) follows from the normalization
$\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\Omega}\sigma=1$ and the definition of
$\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int$.

We now introduce the nondimensional temperature

$u^{*}\equiv{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\dfrac{T^{*}-T_{\text{ref}}}{(\ell
q)/k_{\inf}}}}$

substitution of which into equations (47) – (49) yields

$\displaystyle-\boldsymbol{\nabla}\cdot(\kappa\boldsymbol{\nabla}u^{*})$
$\displaystyle=\sigma\gamma\text{ in }\Omega,$ (52)
$\displaystyle\kappa\partial_{n}u^{*}$ $\displaystyle=-1\text{ on
}\partial\Omega,$ (53) $\displaystyle\int_{\Omega}\sigma u^{*}$
$\displaystyle=0.$ (54)

The average resistance then takes the form

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{R}_{\text{avg}}}}$
$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\equiv\dfrac{\ell}{k_{\inf}|\partial\Omega|}}}\left(-\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\partial\Omega}u^{*}\right),$ (55)

in which we have taken advantage of (54) to eliminate the integral over
$\Omega$. We now choose to write

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{R}_{\text{avg}}=\dfrac{\mathscr{L}_{\text{cond}}}{k_{\inf}|\partial\Omega|}}},$
(56)

which implicitly defines the nondimensional conduction resistance length scale

$\displaystyle\mathscr{L}_{\text{cond}}\equiv{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{L}_{\text{cond}}/\ell}}$
$\displaystyle=-\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\partial\Omega}u^{*}.$ (57)

We now note from inspection of equations (52) – (54) that
$u^{*}=|\Omega|^{1/2}\psi^{\prime 0}$ for $\psi^{\prime 0}$, the solution to
(44) – (46). It then further follows from the weak form of (44) – (46)
(provided in Appendix A.4), definition of $\phi$ (37), and definition of
$\mathscr{L}_{\text{cond}}$ (57) that

$\displaystyle\phi$
$\displaystyle=-|\Omega|^{-1/2}\int_{\partial\Omega}\psi^{\prime 0}$
$\displaystyle=-|\Omega|^{-1}|\partial\Omega|\,\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\partial\Omega}u^{*}$
$\displaystyle=\mathscr{L}_{\text{cond}}/\mathscr{L}\,,$ (58)

and hence we may write

$\displaystyle\mathscr{L}_{\text{cond}}$ $\displaystyle=\phi\,\mathscr{L}.$
(59)

We thus observe that $\mathscr{L}_{\text{cond}}$, the correct conduction
length scale to represent the body average thermal resistance, is the usual
length scale $\mathscr{L}$ amplified by $\phi$. Note that
$\mathscr{L}_{\text{cond}}$ will depend, through $\phi$, on $\sigma$,
$\kappa$, and of course $\Omega$. Our discussion here thus extends, in a
quantitative fashion, the usual notion of resistance between two surfaces to
the case of resistance from body “center” to body boundary. Furthermore, we
note that the definition yields a positive length scale:

######  Proposition 4.1 (Positive Conduction Length Scale,
$\mathscr{L}_{\text{cond}}$).

For any given dunking problem with admissible parameters $\sigma$, $\kappa$,
$B$, and $\Omega$; the corrected conduction length scale

$\displaystyle\mathscr{L}_{\text{cond}}>0.$ (60)

This proposition can be proven by considering (58), positivity of
$\mathscr{L}$ through positive $|\Omega|$ and ${|\partial\Omega|}$, and $\phi$
being a strictly positive quantity — which we shall later demonstrate.

As we shall discuss below, for many “typical” domains (slab, cylinder,
sphere), $\phi\approx\pazocal{O}(1)$; however, for other domains, for example
the finned block shown in Figure 1, $\phi$ and thus
$\mathscr{L}_{\text{cond}}/\mathscr{L}$ may be large. This large value of
$\phi$ reflects a conjoined domain with disparate length scales associated
with each of the component domains — often found in engineering applications.
We investigate the dependence of $\phi$ on geometry and composition in Section
5.

$2\delta$ Figure 1: Fin attached to a square block.

###  4.2 Thermal Circuit

In this section, we continue in order to motivate approximations which will be
rigorously analyzed: the main result being a ‘corrected’ ordinary differential
equation for $u_{\text{avg}}$, leading to the second-order lumped
approximation. To that effect, we now introduce a thermal circuit, shown in
Figure 2, with three nodes and associated nodal temperatures,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{T}_{\text{avg}}}}$,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{T}_{{\partial{\text{avg}}}}}}$,
and
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{T}_{\infty}(=T_{\infty})}}$.
The nodes associated with
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{T}_{\text{avg}}}}$
and
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{T}_{{\partial{\text{avg}}}}}}$
are connected by a body average resistance
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{R}_{\text{avg}}}}$,
and the nodes
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{T}_{{\partial{\text{avg}}}}}}$
and
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{T}_{\infty}}}$
are connected by a standard “heat transfer coefficient” resistance,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{R}_{h}\equiv
1/(h|\partial\Omega|)}}$. It is important to note that $\hat{\cdot}$ refers to
approximation to the standard quantities:
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{T}_{\text{avg}}}}$
and
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{T}_{\partial{\text{avg}}}}}$
are representative of $T_{\text{avg}}$ and $T_{\partial{\text{avg}}}$, but in
general not equal to those quantities.

${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{T}_{\text{avg}}}}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{T}_{\partial{\text{avg}}}}}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{T}_{\infty}}}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{R}_{\text{avg}}}}$${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{R}_{h}}}$
Figure 2: 3-node thermal circuit.

We may now apply standard linear circuit theory to this network. In
particular, we can evaluate the heat transfer rate through our circuit as

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{Q}=\dfrac{\hat{T}_{\text{avg}}-\hat{T}_{\infty}}{\mathscr{R}_{\text{eq}}}}}$
(61)

and the temperature
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{T}_{{\partial{\text{avg}}}}}}$
as

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\dfrac{\hat{T}_{{\partial{\text{avg}}}}-\hat{T}_{\infty}}{\hat{T}_{\text{avg}}-\hat{T}_{\infty}}=\dfrac{\mathscr{R}_{h}}{\mathscr{R}_{\text{eq}}}}}$
(62)

where
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{R}_{\text{eq}}}}$
is the (total, or) equivalent thermal resistance of the circuit,

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{R}_{\text{eq}}\equiv\mathscr{R}_{\text{avg}}+\mathscr{R}_{h}=}}\,\dfrac{\phi{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{L}}}}{{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{k_{\inf}|\partial\Omega|}}}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\,+\dfrac{1}{h|\partial\Omega|}}}.$
(63)

We thus observe that $\phi$ large increases the body average resistance and
thus decreases the heat transfer rate. From another perspective, $\phi$ large
decreases the fraction of the temperature difference
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{(\hat{T}_{\text{avg}}-\hat{T}_{\infty})}}$
associated with
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{R}_{h}}}$
which thus again decreases the heat transfer rate. Note our analysis is
perforce approximate since our average resistance assumes uniform source and
boundary flux; however, we anticipate that our approximation will be quite
good since for small Biot number, the temperature distribution — for example,
our first eigenmode — is quite uniform.

We now note that

$\displaystyle\text{Bi}^{\prime}_{\text{dunk}}\equiv\phi\text{Bi}_{\text{dunk}}={\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\dfrac{\mathscr{R}_{\text{avg}}}{\mathscr{R}_{h}}}}$
(64)

which is precisely the usual definition of the Biot number in terms of
resistance ratio — but now with the improved conduction length scale provided
by (59). We can then express
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{R}_{\text{eq}}}}$
as

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{R}_{\text{eq}}=\dfrac{1}{h|\partial\Omega|}}}(1+\text{Bi}^{\prime}_{\text{dunk}}).$
(65)

Our interest in (65) is in the limit of small Biot number. In the limit of
large Biot number we would write instead
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{R}_{\text{eq}}}}=\phi{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\left[\mathscr{L}/(k_{\inf}|\partial\Omega|)\right]}}\left[1+(\text{Bi}^{\prime}_{\text{dunk}})^{-1}\right]$.

We now turn to the time-dependent problem. We associate to the
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{T}_{\text{avg}}}}$
node a heat capacitance
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{(\overline{\rho
c})|\Omega|}}$ which then directly yields

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{(\overline{\rho
c})|\Omega|\partial_{t}\hat{T}_{\text{avg}}=-\dfrac{\hat{T}_{\text{avg}}-\hat{T}_{\infty}}{\mathscr{R}_{\text{eq}}}}}.$
(66)

We now define
$\hat{u}^{*}_{\text{avg}}\equiv{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{(\hat{T}_{\text{avg}}-\hat{T}_{\infty})/(T_{\text{i}}-\hat{T}_{\infty})}}$
and nondimensionalize (66) to obtain

$\displaystyle\partial_{t}\hat{u}^{*}_{\text{avg}}=-\dfrac{B\gamma}{1+\phi
B/\gamma}\hat{u}^{*},\quad t>0$ (67)

subject to $\hat{u}^{*}_{\text{avg}}(t=0)=1$ from our uniform initial
condition. We can also identify the eigenproblem associated with (67):

$\displaystyle\dfrac{B\gamma}{1+\phi
B/\gamma}\hat{\psi}=\hat{\lambda}\hat{\psi},$ (68)

where $(\hat{\psi},\hat{\lambda})$ is the single eigenpair associated with our
scalar first-order problem.

###  4.3 Lumped Approximations

We begin with the eigenproblem in order to arrive at the lumped approximation.
We directly identify from (68) our first-order approximation — neglecting the
$\pazocal{O}(B^{2})$ contributions — as $\tilde{\lambda}^{1}=B\gamma$, and our
second-order approximation — including now the $\pazocal{O}(B^{2})$ terms — as
$\tilde{\lambda}^{2\text{P}}=\frac{B\gamma}{1+\phi B/\gamma}$. We have already
presented and analyzed these approximations in Proposition 3.4, respectively.
We can now also include a physical interpretation in particular for the
second-order approximation: for larger $\phi$, hence larger conduction length
scale, the heat transfer coefficient resistance is relatively smaller compared
to the conduction (average) resistance, and thus
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\hat{T}_{{\partial{\text{avg}}}}}}$
is closer to
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T_{\infty}}}$;
thus the heat transfer rate to the environment,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{h|\partial\Omega|(T_{\text{avg}}-T_{\infty})}}/(1+\text{Bi}^{\prime}_{\text{dunk}})$
is reduced; and hence the time constant — $1/\tilde{\lambda}^{2\text{P}}$ — is
increased.

Turning now to the initial value problem, we can directly identify from (67)
our first-order approximation for the domain average — neglecting the
$\pazocal{O}(B^{2})$ contributions — as
$\tilde{u}_{\text{avg}}^{1}\equiv\exp(-B\gamma t)$, the ‘classic’
approximation — and then the second-order approximation for the domain average
— now including the $\pazocal{O}(B^{2})$ terms — as
$\tilde{u}_{\text{avg}}^{2}\equiv\exp\left(-\frac{B\gamma t}{1+\phi
B/\gamma}\right)$; note again that larger $\phi$, and larger conduction length
scale, leads to longer time scales for equilibration. Finally, it follows from
(62) that $\tilde{u}_{\Delta}^{2}\equiv\frac{\phi B}{\gamma}$ (independent of
time) — we use superscript $2$ in reference to its derivation from a second-
order approximation. We summarize these three results in Definition 4.2. We
also note that alternative approaches exist to derive the lumped approximation
$\tilde{u}_{\text{avg}}^{1}$; one such approach, Galerkin projection of the
PDE, is presented in Appendix A.6: closely related to the Rayleigh quotient
(Appendix C.4) for a constant function.

######  Definition 4.2 (Lumped Approximations: Summary).

The first and second-order lumped approximations to the QoI, $u_{\text{avg}}$
and $u_{\Delta}$, are given by

$\displaystyle\tilde{u}_{\text{avg}}^{1}(t;B)$
$\displaystyle\equiv\exp(-B\gamma t)\;,$ (69)
$\displaystyle\tilde{u}_{\text{avg}}^{2}(t;B)$
$\displaystyle\equiv\exp\left(-\dfrac{B\gamma t}{1+\phi B/\gamma}\right)\;,$
(70) $\displaystyle\tilde{u}_{\Delta}^{2}(B)$ $\displaystyle\equiv\dfrac{\phi
B}{\gamma}\;.$ (71)

The associated error analysis for these approximations is provided in Section
6.

We note that a second-order lumped approximation is proposed already in [4]
based on the eigenvalue expansion (in [4], for the particular case of a
homogeneous sphere). This solution based on second-order Taylor approximation
of the principal eigenvalue is different from our solution,
$\tilde{u}_{\text{avg}}^{2}$, based on the second-order $[1,1]$ Padé
approximation [10] to the principal eigenvalue: the two proposals agree to
second order in $B$ as $B$ tends to zero; however, the Taylor approximation is
only meaningful for $\frac{\phi B}{\gamma}<1$ whereas the Padé approximation
yields stable results for any $B$. It is interesting that the engineering
resistance arguments directly yield the Padé approximation; see also [11]. We
also note that $\tilde{u}_{\Delta}^{2}$ derived from engineering arguments is
equivalent to the first-order truncation of $u_{\Delta}$ limit for large $t$:

$\displaystyle\lim_{t\rightarrow\infty}u_{\Delta}(t;B)=\frac{M(\psi)-H(\psi)}{M(\psi)}=0+\underbrace{\frac{M(\psi^{\prime\,0})-H(\psi^{\prime\,0})}{M(\psi^{0})}B}_{\tilde{u}_{\Delta}^{2}(B)}+\,\pazocal{O}(B^{2})\,.$
(72)

It is also important that calculation of the dimensional quantity
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{u_{\text{avg}}}}$
requires only ${u}_{\text{avg}}$,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\ell}}$,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\overline{\rho
c}}}$, and
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\overline{h}}}$.
Hence again we need only the average, not local, heat transfer coefficient. It
is the average heat transfer coefficient,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\overline{h}}}$,
and not the local heat transfer coefficient,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{h}}$,
which is most readily measured and most easily archived. Furthermore, for
heterogeneous bodies,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\overline{\rho
c}}}$ can be evaluated solely on the basis of constituent properties and
volume fractions — the details of the material distribution are not required.

##  5 Solutions of the Sensitivity Equation: $\phi$

We now turn to a detailed analysis of $\psi^{\prime\,0}$ and $\phi$. The
quantity $\phi$ appears in our second-order approximations and in all our
error estimates.

We first introduce some notation relevant to our study of the geometric
dependence of $\phi$. Let

$\displaystyle\Xi(x_{\text{pg}}^{1}\in\mathbb{R}^{2},x_{\text{pg}}^{2}\in\mathbb{R}^{2},\ldots,x_{\text{pg}}^{n_{\text{pg}}}\in\mathbb{R}^{2})$
(73)

denote the polygon $\Omega_{\text{pg}}\subset\mathbb{R}^{2}$ with boundary

$\displaystyle\partial\Omega_{\text{pg}}\equiv\bigcup_{i=1}^{n_{\text{pg}}}\;\overline{x_{\text{pg}}^{i}\;x_{\text{pg}}^{i+1}}\;.$
(74)

We shall assume that the vertices $x_{\text{pg}}^{i},1\leq i\leq
n_{\text{pg}}$, are chosen such that the interiors of the segments
$\overline{x_{\text{pg}}^{i}\;x_{\text{pg}}^{i+1}},1\leq i\leq n_{\text{pg}}$,
are mutually non-intersecting. The polygons are simply connected but not
necessarily convex or star-shaped.

###  5.1 General Properties of $\phi$

We summarize here several properties of $\phi(\sigma,\kappa;\Omega)$ which can
be gainfully applied in the evaluation of $\phi$ and the application of $\phi$
to error analysis:

######  Proposition 5.1 ($\phi$ is Strictly Positive).

For any given admissible $\sigma$ and $\kappa$, $\phi(\sigma,\kappa)>0$.

######  Proposition 5.2 ($\phi$ is Scale-Invariant).

We are given $\Omega_{1}\subset\mathbb{R}^{d}$ and a translation-rotation-
uniform dilation map of $\Omega_{1}$,
${\Omega}_{2}\equiv\\{\mathbb{M}(x)\,|\,x\in\Omega_{1}\,\\}$; here
$\mathbb{M}(x)=\mathbb{T}+\mathbb{Q}\,x+\alpha\,\mathbb{I}\,x$, where $\alpha$
is a positive real scalar, $\mathbb{T}\in\mathbb{R}^{d}$ is a translation
vector, $\mathbb{Q}\in\mathbb{R}^{d\times d}$ is a rotation matrix, and
$\mathbb{I}\in\mathbb{R}^{d\times d}$ is the identity. For given
$\sigma\equiv\sigma_{1},\kappa\equiv\kappa_{1}$, $\Omega\equiv\Omega_{1}$ we
evaluate $\phi_{1}\equiv\phi(\Omega_{1})$ from Proposition 3.5; for given
$\sigma\equiv\sigma_{2}\equiv\sigma_{1}\circ\mathbb{M}^{-1},\kappa\equiv\kappa_{2}\equiv\kappa_{1}\circ\mathbb{M}^{-1}$,
$\Omega\equiv\Omega_{2}$ we evaluate $\phi_{2}\equiv\phi(\Omega_{2})$. Then,
for any $\alpha$, $\mathbb{T}$, and $\mathbb{Q}$, $\phi_{1}=\phi_{2}$.

######  Proposition 5.3 ($\phi$ is Monotonic in $\kappa$).

For $\kappa_{1}$ and $\kappa_{2}$ such that $\kappa_{1}\leq\kappa_{2},\forall
x\in\Omega$, $\phi(\sigma,\kappa_{1})\geq\phi(\sigma,\kappa_{2})$.

We henceforth indicate the dependence of $\phi$ as
$\phi(\sigma,\kappa;\Omega_{\text{ref}})$ and, from Proposition 5.2, we can
apply this value of $\phi$ to any domain $\Omega$ which is translation-
rotation-uniform dilation map of $\Omega_{\text{ref}}$. We can further state,
from Proposition 5.3, that
$\phi(\sigma,\kappa;\Omega_{\text{ref}})\leq\phi(\sigma,1;\Omega_{\text{ref}})$
for any $\kappa\geq 1$; given the latter bound, we shall often restrict
attention to $\kappa=1$ even if the result can be developed for general
$\kappa$. We provide in Appendix C.8 sketches of the proofs of Propositions
5.1 – 5.3.

Proposition 5.3 will permit us to develop our error estimators even for the
case — quite common in engineering practice — of heterogeneous material
composition. In the next proposition we provide also some insight into the
dependence of $\phi$ on $\sigma$ and in particular on deviations of $\sigma$
from uniform:

######  Proposition 5.4 ($\phi$ Dependence on $\sigma$).

We consider here the case in which $\kappa=1$ and we furthermore assume that
the domain $\Omega$ is fixed; we may thus abbreviate $\phi$ as simply
$\phi(\sigma)$. Let us introduce

$\displaystyle\mu\equiv\inf_{w\in Z_{0}(1)}\dfrac{a_{0}(w,w;1)}{m(w,w;1)}\;.$
(75)

We further define

$\displaystyle\delta(\sigma)\equiv\left(\dfrac{\gamma^{2}}{\mu}\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\Omega}(\sigma-1)^{2}\right)^{1/2}\;,$ (76)

in terms of which we may then bound $\phi(\sigma)$ as

$\displaystyle\phi(\sigma)\leq\phi_{\text{UB}}(\sigma)\equiv\left(\left(\phi(1)\right)^{1/2}+\delta(\sigma)\right)^{2}\;,$
(77)

where $\phi(1)$ refers to $\phi$ for $\sigma=1$. We recall that
$\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int$ refers to the average of the integrand over the
domain.

We provide in Appendix C.9 a proof of Proposition 5.4.

We note that $\gamma^{2}/\mu$ is scale invariant. We further note that $\mu$
is simply the second (or any second) eigenvalue of the Neumann Laplacian
eigenproblem on the domain $\Omega$. It thus follows that $\mu$ may be
estimated very efficiently, and often with sufficient precision, by
application of the Payne-Weinberger [12] lower bound for convex domains

$\displaystyle\mu\geq\mu_{\text{LB}}\equiv\pi^{2}\mathscr{D}_{\Omega}^{-2}\;$
(78)

where

$\displaystyle\mathscr{D}_{\Omega}=\sup_{x\in\Omega,y\in\Omega}|x-y|\;$ (79)

requires calculation of only the diameter of our domain. For later numerical
results in Section 5.4, we shall numerically evaluate $\mu_{h}(=\mu)$; in
general, $(\cdot)_{h}$ shall refer to finite element approximation of
$(\cdot)$. In actual practice, the lower bound $\mu_{\text{LB}}$ is more
convenient.

Proposition 5.4 demonstrates that in fact there are many physical situations
in which $\phi(\kappa=1,\sigma=1)$ — easily determined without detailed
knowledge of the heterogeneous material composition — can serve as a surrogate
for $\phi(\kappa=1,\sigma)$. We note in particular that either if the two
materials are similar in thermophysical properties or if the two materials
occupy commensurate volume (fraction), then
$\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\Omega}(\sigma-1)^{2}$ is either small or of order
unity. As an example of the latter, we consider two materials of very
different properties each of which occupies half of $\Omega$:

######  Proposition 5.5 (Equi-Volume-Fraction Two-Phase Composite).

We consider two materials (of different thermophysical properties) each of
which constitutes one half of the volume of our domain $\Omega$. Then
$\mathchoice{{\vbox{\hbox{$\textstyle-$
}}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$
}}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$
}}\kern-3.82503pt}}\\!\int_{\Omega}(\sigma-1)^{2}\leq 1$.

Note that even in this equi-volume-fraction case, $\delta(\sigma)$ may still
be large due to geometry effects through $\mu$. We provide in Appendix C.10 a
proof of Proposition 5.5.

###  5.2 Closed-Form Results

We first present in Table 1 the values of
$\phi(\kappa=1,\sigma=1;\Omega_{\text{ref}})$ for the canonical reference
domains: interval in $\mathbb{R}^{1}$, disk in $\mathbb{R}^{2}$, and sphere in
$\mathbb{R}^{3}$. We choose $\kappa=1$ and $\sigma=1$ for two reasons: first,
often in practice we encounter homogeneous properties; second, thanks to
Proposition 5.3 and Proposition 5.4, we can bound the effect of deviations
from uniform properties in terms of readily measured quantities. Note the
sphere calculation is presented in [4] for the case of homogeneous properties
($\sigma=\kappa=1$) by direct eigenvalue perturbation expansion; in Appendix B
we reproduce this homogeneous sphere result based on the (strong form of the)
Sensitivity Equation (Proposition 3.5).

${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\Omega_{\text{ref}}\equiv\\{x\in\mathbb{R}^{d}\,|\,|x|\leq\ell\\}}}$ | $\phi$
---|---
$\quad\quad d=1$: Interval | 1/3
$\quad\quad d=2$: Disk | 1/2
$\quad\quad d=3$: Sphere | 3/5
Table 1: Values of $\phi(\kappa=1,\sigma=1;\Omega_{\text{ref}})$ for Canonical
Domains in Dimension $d\in\\{1,2,3\\}$. Note that we include a nominal length
scale for each geometry, however it follows from scale invariance, Proposition
5.1, that $\phi(\kappa=1,\sigma=1;\Omega_{\text{ref}})$ is independent of
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\ell}}$
(Proposition 5.2).

It is also possible to analyze selected triangles in closed form. We present
here a limiting case for a right triangle which is algebraically simple but
illustrates the key point:

######  Proposition 5.6 (Evaluation of $\phi$ for Selected Triangles).

We consider the right triangle $\Omega^{\text{right
triangle}}\subset\mathbb{R}^{2}$ given by

$\displaystyle\Xi^{\text{right
triangle}}(W)(x_{\text{pg}}^{1}\equiv(0,0),x_{\text{pg}}^{2}\equiv(W,0),x_{\text{pg}}^{3}\equiv(0,1))\;,$
(80)

for $W\in\mathbb{R}_{+}$. Then

$\displaystyle\phi(\kappa=1,\sigma=1;\Xi^{\text{right
triangle}}(W))\sim\dfrac{2}{3}W^{-2}\text{ as }W\rightarrow 0\;.$ (81)

As by-products of this analysis, we find that

$\displaystyle\phi(\kappa=1,\sigma=1;\Xi^{\text{right
triangle}}(1))=\frac{4}{3}$ (82)

and

$\displaystyle\phi(\kappa=1,\sigma=1;\Xi^{\text{equilateral triangle}})=1\;.$
(83)

Recall that, as always, these results are scale invariant (Proposition 5.2).

We provide a proof of Proposition 5.6 in Appendix C.11.

We note that, whereas the canonical domains of Table 1 all yield $\phi$ on the
order of unity, a triangular domain can, in the limit $W\rightarrow 0$,
produce a very large $\phi$. In the next section we will encounter other
domains (and thermophysical property limits) which exhibit large $\phi$. And
we subsequently demonstrate, in the context of error analysis, that the
standard engineering criterion for applicability of the lumped approximation —
$\text{Bi}_{\text{dunk}}$ small — implicitly assumes
$\phi\approx\pazocal{O}(1)$. This relates to the physical (thermal circuit)
analysis of Section 4.2 (i.e.,
$\phi\approx\pazocal{O}(1)\rightarrow\text{Bi}_{\text{dunk}}\approx\text{Bi}^{\prime}_{\text{dunk}}$).

Finally, $\phi$ may be readily evaluated for tensorized domains. We provide
here a particular but illustrative case (and for simplicity, we consider
$\kappa=\sigma=1)$:

######  Proposition 5.7 ($\phi$ for Tensorized Domains).

We introduce, for a given two-dimensional domain
$\Omega^{\text{2D}}\subset\mathbb{R}^{2}$, the tensorized (or extruded) three-
dimensional domain
$\Omega\subset\mathbb{R}^{3}\equiv\Omega^{\text{2D}}\times(0,L)$. Then

$\displaystyle\phi(1,1;\Omega)=\phi_{\text{2D}}(1,1;\Omega^{\text{2D}})+\dfrac{1}{3}\;,$
(84)

where $\phi_{\text{2D}}$ refers to the value of $\phi$ for a two-dimensional
domain (in our case here, $\Omega^{\text{2D}}$). As always, our result is
scale invariant with respect to $\Omega^{\text{2D}}$ (Proposition 5.2).

We provide a proof of Proposition 5.7 in Appendix C.12. Note that Proposition
5.7 may be extended to non-uniform properties — and in particular we may
appeal to Proposition 5.4 — if we require
$\sigma(x)=\sigma^{\text{2D}}(x_{1},x_{2})\sigma^{\text{1D}}(x_{3})$ for
$\sigma^{\text{2D}}$ and $\sigma^{\text{1D}}$ unity mean. Note that $\phi$ of
Proposition 5.7 does not depend on $L$, and in particular there is no
distinction between small $L$ and large $L$. This result is counter-intuitive,
however it is important to recall that $\phi$ is a coefficient of the second-
order contribution to the eigenvalue. Recall also our assumption that the heat
transfer coefficient, Robin coefficient, is uniform over the entire boundary
$\partial\Omega$.

We consider two illustrative applications of Proposition 5.7 in concert with
respectively Table 1 and Proposition 5.6: for a cylinder of radius $R$ and
length $L$, $\phi=1/2+1/3=5/6$; for a prismatic right triangle
$\Xi^{\text{right triangle}}(W)\times(0,L)$, $\phi\sim\frac{2}{3}W^{-2}$ as
$W\rightarrow 0$. A small variation of Proposition 5.7 (in conjunction with
Table 1) readily reveals that, for a rectangle of dimensions $L_{1}\times
L_{2}$ and an orthogonal parallelepiped of dimensions $L_{1}\times L_{2}\times
L_{3}$, $\phi=2/3$ and $\phi=1$, respectively.

###  5.3 Foundations for Alternative Evaluation Strategies

As discussed in the introduction, our lumped approximation, a simple
exponential, has many attractive features: transparency — explicit parameter
dependence; flexibility — ready incorporation into design or optimization or
system analysis; rapid evaluation — just a few standard elementary function
evaluations.

However, our error estimates shall require the new quantity $\phi$: we already
observe in Proposition 3.3 that the relative error estimate for the first-
order lumped approximation to the first eigenvalue is given by
$\text{Bi}^{\prime}_{\text{dunk}}(\equiv\phi\,\text{Bi}_{\text{dunk}})$; we
shall similarly observe in the Section 6 that the error in the first-order
lumped approximation to $u_{\mathrm{avg}}$ is similarly given by
$\text{Bi}^{\prime}_{\text{dunk}}/e$ (Proposition 6.5). We now understand that
$\phi$ may, in fact, be large — and hence the implicit engineering assumption
(say, $\phi=1$) can not, in general, be trusted: $\phi$ must be evaluated. We
further know that, in general, $\phi$ cannot be obtained in closed form except
for a few simple geometries. We conclude that we must, explicitly or
implicitly, resort to computational treatment of the partial differential
equation (129) for $\psi^{\prime 0}$ in order to evaluate $\phi$ from (130).
Do we thus lose the advantages of the lumped approximation? In the next
paragraphs we show that, in fact, we do not lose these advantages.

We first note that although the set of geometries amenable to closed-form
evaluation of $\phi$ (Table 1, Proposition 5.6) is small, these geometries are
of practical interest as they often arise in both engineered and natural
systems, and furthermore, these geometries can serve to approximate more
general shapes. We also emphasize that this set of geometries is further
enriched by Proposition 5.7, providing closed-form $\phi$ expressions for
finite circular cylinder, orthogonal parallelepiped, and certain triangular
prisms.

However, complex geometries not amenable to simple approximations often arise
in engineering applications. In practice, we do not lose the advantages of the
lumped approximation: a sufficiently accurate approximation to our single
elliptic partial differential equation can be very rapidly obtained by
adaptive finite element methods (see Section 5.4); furthermore, more recent
software advances, in particular containers and microservices, permit cloud-
based deployment to virtually any user — with no issues of licenses or
installation or compilation or access. Nevertheless, the lumped approximation
requires only the volume and surface area of the body, whereas our partial
differential equation requires a full geometry (and material composition)
description — often an inconvenience. We may thus wish to supplement direct
evaluation of $\phi$ with a more convenient alternative strategy in which
$\phi$ is replaced by an inexpensive surrogate.

As regards geometry, there are two options. In both cases, we are given the
desired domain, $\Omega$, and a dictionary of $(\Omega,\phi)^{k}$ pairs. In
the first option, we estimate $\phi$ as $\phi^{k}$ for the element $k$ of the
dictionary for which $\text{dist}(\Omega,\Omega^{k})$ (with suitable
registration) is smallest. In the second option, we identify a few features
which permit approximate separation of our dictionary entries for $\phi$
smaller and $\phi$ larger; in the offline stage we train the classifier, and
then in the online stage we query the classifier for any given desired domain,
$\Omega$. We are thus presented with two mathematical questions: (for the
first option) How should we choose $\mathrm{dist}$, in particular to provide
Lipschitz continuity for $\phi$ — $\mathrm{dist}(\Omega^{i},\Omega^{j})$ small
implies $|\phi^{i}-\phi^{j}|$ small [13]?; (for the second option) How should
we choose our features to provide adequate separation between the $\phi$
smaller and $\phi$ larger subsets of the dictionary (or, more generally,
training set)?

We cannot yet provide rigorous answers. However, we have conjectures which we
will then support (in the next section) with computational evidence. A
subsequent manuscript will elaborate on the material presented here and
provide (and assess) specific surrogate techniques.

  * For the distance we propose

$\displaystyle\text{dist}(\Omega^{i},\Omega^{j})\equiv
C_{1}d_{\mathrm{Hausdorff}}(\Omega^{i},\Omega^{j})+C_{2}\left||\partial\Omega^{i}|-|\partial\Omega^{j}|\right|,$
(85)

where $d_{\mathrm{Hausdorff}}$ is the Hausdorff distance [14] and $C_{1}$ and
$C_{2}$ are positive real constants.

  * The geometric features should, in general, increase with $\phi$. For the first feature, we propose

$\displaystyle\text{Feature}_{1}$ $\displaystyle\equiv\text{Modified Weighted
Inverse Angle (MWIA)},$ (86)

$\text{Feature}_{1}$ is complicated and we thus restrict attention to
$\Omega\equiv\Omega^{\text{2D}}\subset\mathbb{R}^{2}$ (or associated extruded
domains $\Omega\times(0,W)$ in $\mathbb{R}^{3}$) and furthermore polygonal
domains $\Omega^{\text{2D}}$. In $\text{Feature}_{1}$ we consider all
polygonal domains which are sufficiently close to $\Omega^{\text{2D}}$ in
dist; for each such sufficiently close domain, we sum the inverse of the
squared angles over all vertices weighted by the relative perimeter of the
edges which share the vertex. The ‘modified’ qualifier indicates that a pre-
processing procedure in which small segments (compared to its neighbors) are
removed to account for ‘smoothed’ sharp angles which may fool an ‘unmodified’
procedure. In short, $\text{Feature}_{1}$ will be large if some good fraction
of the boundary corresponds to small interior angles.

  * For the second feature, we propose

$\displaystyle\text{Feature}_{2}$ $\displaystyle\equiv\text{InRadius
$\gamma$}.$ (87)

$\text{Feature}_{2}$ is much simpler than $\text{Feature}_{1}$, since the
InRadius, the radius of the largest inscribed ball, is readily evaluated or
estimated. In short, $\text{Feature}_{2}$ will be large if there are extended
regions which contribute substantially to $|\partial\Omega|$ but not to
$|\Omega|$.

Examples will be provided in the next section.

The choice of distance dist and also $\text{Feature}_{1}$ and
$\text{Feature}_{2}$ can be understood in terms of the simple resistance
arguments presented in Section 4. As an example, consider two domains
illustrated in Figure 3: $\Omega^{1}$ is the unit-radius disk; $\Omega^{2}$ is
the unit-radius disk but with a highly oscillatory (in space) boundary — a
version of the latter, a gear, shall be studied in the next section. We
observe that for $\Omega^{2}$, $\mathscr{L}=1/\gamma$ will be very small, even
though clearly the relevant conduction length scale is the radius (= 1): to
reflect the latter, $\phi$ for $\Omega^{2}$ must be large compared to
$\phi=1/2$ for our disk, $\Omega^{1}$. And, as we would hope,
$\text{dist}(\Omega^{1},\Omega^{2})$ will be large — thanks to the term
associated with the difference in boundary measure. And, again as we would
hope, for $\Omega^{1}$, both $\text{Feature}_{1}$ and $\text{Feature}_{2}$ are
not large, where as for $\Omega^{2}$, both $\text{Feature}_{1}$ and
$\text{Feature}_{2}$ are large. We shall observe other examples in which only
one of the two features will be large — hence the need for both features.

Figure 3: Illustrative domains $\Omega^{1}$ and $\Omega^{2}$.

To address material heterogeneity — non-uniform thermophysical properties —
our Proposition 5.4 suggests very simple bounds that in practice can yield
sufficiently sharp results; we shall provide numerical examples in the next
section demonstrating the effectiveness of the proposed bounds.

###  5.4 Representative Computational Results

In this section we provide some illustrative computational results. We exploit
standard $\mathbb{P}_{2}$ finite element approximation with adaptive
refinement [7] and extrapolation based error estimation. In the case of curved
geometries we consider isoparametric approximation and uniform refinement. In
the below, we report the finite element approximation $\phi_{h}$ and also the
associated finite element error estimator $e^{\phi}_{h}$; in practice (and in
theory as $h\rightarrow 0$), $|\phi-\phi_{h}|\lessapprox e^{\phi}_{h}$. The
error estimator is based on the solution convergence rate for a successively
mesh refinement [15, 16]. In all cases we consider domains
$\Omega_{\mathrm{2D}}\subset\mathbb{R}^{2}$, but recall we can readily extend
these results to extruded three-dimensional domains thanks to Proposition 5.7.
Note the finite element approximation is based on the saddle formulation of
Appendix A.5: the latter effectively incorporates the constraint associated
with $Z_{0}(\sigma)$ and furthermore admits a sparse solution.

We shall first focus on the dependence of $\phi$ on $\Omega$; we set
$\sigma=\kappa=1$. We summarize our results, in Table 2 and associated Figures
4 – 11, for a variety of geometric configurations: labels in Table 2 will be
used to reference test domains and the acronyms will be expanded as a
pneumonic in associated figures. We observe, consistent with the asymptotic
theory of Proposition 5.6, that $\phi$ for SART (Figure 4) is quite large (and
indeed very close to the asymptotic result). We also observe in Figure 4 the
rapid decrease of $\psi^{\prime 0}$ and associated large gradients of
$\psi^{\prime 0}$ into the corner; the latter gives rise to the large value of
$\phi$. But we then can conclude from SARTF (Figure 5) that the sharp angle is
not required for large $\phi$ (less than 1% $\phi$ difference compared to
SART); note the fillet relieves the sharp corner but only modifies an
extremely small portion of the domain. We further conclude from SARTC (Figure
6) that neither an acute angle nor high curvature is required for large
$\phi$. We also observe that RECT (Figure 7) yields the analytical $\phi$ for
a rectangle ($2/3$) with numerical error on the order of machine-precision,
and furthermore that a small SART perturbation to RECT, RECTSART (Figure 8),
does not lead to a significant increase in $\phi$. Finally we note that the
two bracket cases, BETA (Figure 9) and BDTA (Figure 10), superficially
somewhat similar, lead to very different values of $\phi$; demonstrating that
a composite domain of small-$\phi$ subdomains can be associated with a large
value of $\phi$. We also remark that in most of the cases, $\phi_{h}/\gamma$
is on the order of $10^{-1}$, meaning $\phi$ and $\gamma$ are correlated;
however, we observe from configuration BDTA that $\phi_{h}$ can be
disproportionately large compared to $\gamma$.

Table 2: Numerical Values of $\phi$ and Related Quantities for Select Domains. Shape | $\phi_{h}$ | $\phi_{h}/\gamma$ | $e_{h}^{\phi}$
---|---|---|---
SART | $9.13\;\\!\cdot\;\\!10^{0}$ | $5.01\;\\!\cdot\;\\!10^{-1}$ | $1.64\;\\!\cdot\;\\!10^{-14}$
SARTF | $9.07\;\\!\cdot\;\\!10^{0}$ | $4.97\;\\!\cdot\;\\!10^{-1}$ | $1.56\;\\!\cdot\;\\!10^{-7}$
SARTC | $9.06\;\\!\cdot\;\\!10^{0}$ | $4.97\;\\!\cdot\;\\!10^{-1}$ | $2.17\;\\!\cdot\;\\!10^{-7}$
RECT | $6.67\;\\!\cdot\;\\!10^{-1}$ | $6.65\;\\!\cdot\;\\!10^{-2}$ | $5.90\;\\!\cdot\;\\!10^{-16}$
RECTSART | $7.04\;\\!\cdot\;\\!10^{-1}$ | $6.98\;\\!\cdot\;\\!10^{-2}$ | $2.63\;\\!\cdot\;\\!10^{-7}$
BETA | $6.95\;\\!\cdot\;\\!10^{-1}$ | $7.60\;\\!\cdot\;\\!10^{-2}$ | $7.45\;\\!\cdot\;\\!10^{-7}$
BDTA | $7.55\;\\!\cdot\;\\!10^{1}$ | $5.07\;\\!\cdot\;\\!10^{0}$ | $8.90\;\\!\cdot\;\\!10^{-4}$

These results which may appear to defy simple summary are, however, consistent
with our conjecture as regards dist (85). In fact, for the cases of Table 2,
${d_{\text{Hausdorff}}}$ suffices for distance. The results of Table 2 are
also consistent with our two features. $\text{Feature}_{1}(\equiv\text{MWIA})$
suffices for all the triangle-related cases; for example it distinguishes
between SART and RECTSART. However, for these triangle-related cases,
$\text{Feature}_{2}(\equiv\text{InRadius }\gamma)$ is not helpful. In
contrast, for discrimination between BETA and BDTA, $\text{Feature}_{1}$ is
not relevant and instead $\text{Feature}_{2}$ is required. We can interpret
the latter in terms of the thermal resistance concept (Section 4): for BDTA,
$\mathscr{L}$ is not representative of the length scale associated with the
conduction thermal resistance $\mathscr{L}_{\text{cond}}$, and must be
corrected by $\phi$; in temporal terms, the small-thickness arm will
equilibrate quickly to ambient and play no further role in the heat transfer
processes.

Figure 4: Geometry Small-Angle Right Triangle, denoted in short as SART:
solution $\psi^{\prime 0}$ and final (adaptively refined) finite element mesh.
The geometry is given by (80) with $W=1/4$. Note we may construct from
$\Omega^{2\mathrm{D}}$ associated extruded domains in $\mathbb{R}^{3}$. Figure
5: Geometry Small-Angle Right Triangle with Fillet (SART with a $10^{-4}$
radius fillet), denoted in short as SARTF: solution $\psi^{\prime 0}$ and
final (adaptively refined) finite element mesh. Note we may construct from
$\Omega^{2\mathrm{D}}$ associated extruded domains in $\mathbb{R}^{3}$. Figure
6: Geometry Small-Angle Right Triangle with Cut ($2\times 10^{-4}$ horizontal
cut modification to SART), denoted in short as SARTC: solution $\psi^{\prime
0}$ and final (adaptively refined) finite element mesh. Note we may construct
from $\Omega^{2\mathrm{D}}$ associated extruded domains in $\mathbb{R}^{3}$.
Figure 7: Geometry Rectangle (width $1/4$ and height $99/100$), denoted in
short as RECT: solution $\psi^{\prime 0}$ and final (adaptively refined)
finite element mesh. Note we may construct from $\Omega^{2\mathrm{D}}$
associated extruded domains in $\mathbb{R}^{3}$. Figure 8: Geometry RECT with
a small SART perturbation (union of RECT and SART), denoted in short as
RECTSART: solution $\psi^{\prime 0}$ and final (adaptively refined) finite
element mesh. Note we may construct from $\Omega^{2\mathrm{D}}$ associated
extruded domains in $\mathbb{R}^{3}$. Figure 9: Geometry Bracket Equi-
Thickness Arms, denoted in short as BETA: solution $\psi^{\prime 0}$ and final
(adaptively refined) finite element mesh. Note we may construct from
$\Omega^{2\mathrm{D}}$ associated extruded domains in $\mathbb{R}^{3}$. Figure
10: Geometry Bracket with Disparate Thickness Arms (ratio 10), denoted in
short as BDTA: solution $\psi^{\prime 0}$ and final (adaptively refined)
finite element mesh. Note we may construct from $\Omega^{2\mathrm{D}}$
associated extruded domains in $\mathbb{R}^{3}$.

We now consider the gear shown in Figure 11(a) for a more comprehensive review
of the effect of $\phi$ as well as further justification of dist (85); note
our extrinsic length scale
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\ell}}$ is
the radius. We consider $m_{\text{teeth}}\geq 4$, even, and define the number
of half teeth as $n=2m_{\text{teeth}}$. We then consider a sequence of
geometries indexed by the number of half teeth, $n$:
$\delta(n;q)\equiv\left(\theta(n)\right)^{q}$, $\theta(n)\equiv\pi/n$; here
$\delta$ and $\theta$ determine the halftooth dimensions in Figure 11(b).
222Our gear example is not realistic for $q<1$ and $n\rightarrow\infty$, not
only in terms of the geometry, but also in the assumption that $B$ is
independent of the number of teeth: in practice, convection and radiation heat
transfer coefficients would decrease as the exterior tooth angle closes. Our
purpose is only to illustrate that our distance dist and features are relevant
to very rough geometries such as the gear for finite $n$. Note that by
symmetry arguments, $\phi$ for the full gear is given by $\phi$ for the half
tooth with homogeneous Neumann conditions on $\Gamma_{2}$ and $\Gamma_{3}$
(Figure 11(a)). The formulation for $\psi^{\prime 0}$ and $\phi$ on the full
gear requires only minor modification to treat the half tooth and associated
homogeneous Neumann conditions: the same geometric factors in equation (129)
yield solvable systems either on the full gear with Robin conditions or on the
half tooth with mixed Robin-Neumann conditions. Of course the interest in the
half tooth is computational expediency.

We present tabulated computational results in Table 3 for the cases of
$(q,n)\in\\{1.6,0.4\\}\times\\{32,256,2048\\}$. For $q=1.6>1$, our gear
approaches a disk (of unit radius) in ${d_{\text{Hausdorff}}}$, in domain
measure $|\Omega^{\mathrm{2D}}|$ and in boundary measure
$|\partial\Omega^{\mathrm{2D}}|$; we observe that $\phi$ for the gear
approaches $\phi$ for the disk while $\text{dist}(gear,disk)$ approaches 0.
For $q=0.4<1$, our gear approaches a disk in ${d_{\text{Hausdorff}}}$ and also
in domain measure $|\Omega^{\mathrm{2D}}|$ but not in boundary measure
$|\partial\Omega^{\mathrm{2D}}|$ (indicated by $\gamma$ progression) — and
hence not in dist; we observe that $\phi$ does not converge to $\phi$ for a
disk. As anticipated, a large change in $\phi$ for two domains is reflected in
a large distance dist between the two domains. Note also that
${d_{\text{Hausdorff}}}$ alone is not sufficient: $\phi$ diverges even as
${d_{\text{Hausdorff}}}$ between successive gears tends to zero. We emphasize
that as we increase the number of teeth, the geometry changes but also the
geometric factors in the equation for $\phi$ (129) change.333We also note that
$\phi/\gamma$, and hence our error bound $\text{Bi}^{\prime}_{\text{dunk}}$,
slowly diverges — and hence so will our bounds for, say, the relative error in
the eigenvalue. Indeed, it can be shown that, for $q<1$, $\lambda_{1}$
converges to a value independent of $B$, observed computationally to be the
Dirichlet eigenvalue for the disk. The gear can also provide insight into our
features. We consider here the case $q<1$. As expected from our discussion of
Section 5.3, both features (86) – (87) diverge as $n$ increases: the gears for
larger $n$ are clearly separated from members of a dictionary which correspond
to small $\phi$.

(a) Full domain and half-tooth domain with boundaries.

(b) Half-tooth computational domain, vertices indicated with $\theta$ and
$\delta$ .

Figure 11: Geometry Gear: full domain and extracted half-tooth ($q=1.6$, $n=32$). We present for the half-tooth the solution $\psi^{\prime 0}_{h}$ and final (adaptively refined) finite element mesh. | $q=1.6$ | $q=0.4$
---|---|---
$n$ | $32$ | $256$ | $2048$ | $32$ | $256$ | $2048$
$\gamma$ | $2.53\;\\!\cdot\;\\!10^{0}$ | $2.05\;\\!\cdot\;\\!10^{0}$ | $2.00\;\\!\cdot\;\\!10^{0}$ | $1.49\;\\!\cdot\;\\!10^{1}$ | $3.91\;\\!\cdot\;\\!10^{1}$ | $1.30\;\\!\cdot\;\\!10^{2}$
$\phi_{h}$ | $7.94\;\\!\cdot\;\\!10^{-1}$ | $5.23\;\\!\cdot\;\\!10^{-1}$ | $5.02\;\\!\cdot\;\\!10^{-1}$ | $3.00\;\\!\cdot\;\\!10^{1}$ | $1.94\;\\!\cdot\;\\!10^{2}$ | $2.12\;\\!\cdot\;\\!10^{3}$
$e^{\phi}_{h}$ | $1.22\;\\!\cdot\;\\!10^{-14}$ | $4.30\;\\!\cdot\;\\!10^{-13}$ | $2.92\;\\!\cdot\;\\!10^{-11}$ | $1.36\;\\!\cdot\;\\!10^{-11}$ | $1.84\;\\!\cdot\;\\!10^{-10}$ | $9.66\;\\!\cdot\;\\!10^{-8}$
Table 3: Numerical results for half-tooth geometries
$(q,n)\in\\{1.6,0.4\\}\times\\{32,256,2048\\}$.

We now turn to a few examples in which the thermophysical properties are non-
uniform, and in particular $\sigma\neq 1$ (we continue to take $\kappa=1$).
The relevant theoretical results are Proposition 5.4 and Proposition 5.5. In
this section we consider a domain $\mathbb{R}^{2}$ ($d=2$) comprised of two
materials, $\sigma_{a}$ over $\Omega_{a}$ and $\sigma_{b}$ over $\Omega_{b}$,
with $\sigma_{a}/\sigma_{b}=1000$ and $|\Omega_{a}|/|\Omega_{b}|=r$ (to be
varied). We consider four configurations: a horizontal interface rectangle
RECTHI (Figure 12), an equal volume-fraction concentric squares EVFCS (Figure
13), and two disparate volume-fraction concentric squares DVFCSLF and DVFCSHF
(Figures 14 – 15) corresponding to light and heavy ‘film’ configurations,
respectively. We present the results for all 4 configurations in Table 4.
Overall, we observe that the $\sigma$ contrast has a significant impact on
$\phi$, and that the bounds of Proposition 5.4 and (78) are effective in some
cases but not in others. We observe that volume fraction has a significant
impact on $\phi$ as demonstrated by the two concentric square examples EVFCS
and DVFCSLF. (Note we do not consider small $\sigma$ contrast since the effect
of non-uniformity, as predicted by Proposition 5.4, is quite small.) We also
note that for the two $r=1$ cases RECTHI and EVFCS, Proposition 5.5 is
applicable and we do indeed observe less than unity $\sigma$ variance (5th
column of Table 4).

Shape | $r$ | $\mu_{h}$ | $\mathchoice{{\vbox{\hbox{$\textstyle-$ }}\kern-7.05008pt}}{{\vbox{\hbox{$\scriptstyle-$ }}\kern-5.31004pt}}{{\vbox{\hbox{$\scriptscriptstyle-$ }}\kern-4.27504pt}}{{\vbox{\hbox{$\scriptscriptstyle-$ }}\kern-3.82503pt}}\\!\int_{\Omega}\left(\sigma-1\right)^{2}$ | $\phi_{\mathrm{UB}}$ | $\phi_{h}(\sigma)$ | $e_{h}^{\phi}$
---|---|---|---|---|---|---
RECTHI | 1 | 9.87 | 0.996 | 15.9 | 8.97 | $1.66\;\\!\cdot\;\\!10^{-14}$
EVFCS | 1 | 9.87 | 0.996 | 4.36 | 1.58 | $7.17\;\\!\cdot\;\\!10^{-6}$
DVFCSLF | 20 | 9.87 | $0.0499$ | 1.21 | 0.732 | $1.77\;\\!\cdot\;\\!10^{-9}$
DVFCSHF | 1/20 | 9.87 | $19.2$ | 40.9 | 0.0181 | $6.07\;\\!\cdot\;\\!10^{-10}$
Table 4: Two-material $\phi_{h}$ results with bounds and associated
quantities. Note: all geometries correspond to $\mu=\pi^{2}$.

We start with a rectangular domain (almost) equivalent to that of RECT shown
in Figure 12. The domain is split into two equal parts along a horizontal
interface ($r=1$). The upper subdomain $\Omega_{a}$ corresponds to
$\sigma_{a}=1.998$ and the lower subdomain $\Omega_{b}$ corresponds to
$\sigma_{b}=1.998\times 10^{-3}$. We observe that the bound for $\phi$ is
within a factor 2 of $\phi_{h}$. We can provide a physical interpretation of
this high-$\phi$ configuration: the lumped approximation assumes uniform decay
to equilibrium over the entire domain; however, it is clear that, for this
configuration, the two subdomains are coupled only along a short interface –
hence effectively decoupled – with disparate $\sigma$-based time scales.
Therefore, a large correction to $\text{Bi}_{\text{dunk}}$ is required.

Figure 12: Geometry Rectangle with Horizontal Interface, denoted in short as
RECTHI: solution $\psi^{\prime 0}$ and final (adaptively refined) finite
element mesh. The geometry is that of RECT with increased height of 1 and
material interface at $y=1/2$. Note we may construct from
$\Omega^{2\mathrm{D}}$ associated extruded domains in $\mathbb{R}^{3}$.

Next we consider EVFCS (shown in Figure 13) for which the two subdomains, an
inner-square ($\Omega_{a}$) and border ($\Omega_{b}$), each have an area of
$1/2$ ($r=1$). The outer region $\Omega_{b}$ corresponds to
$\sigma_{b}=1.998\times 10^{-3}$ and the inner region $\Omega_{a}$ corresponds
to $\sigma_{a}=1.998$. We observe a decrease in the sharpness of the bound for
$\phi$ compared to RECTHI. More importantly, we observe that $\phi_{h}$ is
$\pazocal{O}(1)$, as might be expected on physical grounds: the two subdomains
interact in a coupled fashion, since the thermal flux from $\Omega_{a}$ must
pass through $\Omega_{b}$ prior to leaving the domain. Therefore, compared to
RECTHI, the $\text{Bi}_{\text{dunk}}$ correction is less severe.

Figure 13: Geometry Equi Volume-Fraction Concentric Squares, denoted in short
as EVFCS: solution $\psi^{\prime 0}$ and final (adaptively refined) finite
element mesh. Note we may construct from $\Omega^{2\mathrm{D}}$ associated
extruded domains in $\mathbb{R}^{3}$.

For DVFCSLF (shown in Figure 14), we start with the same domain and topology
as EVFCS, but we increase $r$ to 20, for the same $\sigma$ ratio. Due to this
disparity in the volume-fraction, the $\sigma$ values are significantly
altered to $\sigma_{a}=1.05$ and $\sigma_{b}=1.05\times 10^{-3}$. We do see an
increase in sharpness of the bound for $\phi$ in this case, a result of
dramatically reduced $\sigma$ variance. We also note that, for this
configuration, we expect to see $\phi(\sigma)$ approach that of $\phi(1)=2/3$
as $r$ increases: $\sigma$ variance diminishes while $\gamma$ and $\mu$ are
fixed. In effect, since

$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\frac{(\rho
c)_{a}|\Omega_{a}|}{(\rho
c)_{b}|\Omega_{b}|}}}=\frac{\sigma_{a}|\Omega_{a}|}{\sigma_{b}|\Omega_{b}|}=1000\cdot
20,$ (88)

so material $b$ is effectively absent.

Figure 14: Geometry Disparate Volume-Fraction Concentric Squares with Light
Film, denoted in short as DVFCSLF: solution $\psi^{\prime 0}$ and final
(adaptively refined) finite element mesh. Note we may construct from
$\Omega^{2\mathrm{D}}$ associated extruded domains in $\mathbb{R}^{3}$.

For DVFCSHF (shown in Figure 15), we start with the same domain and topology
as DVFCSLF, but we exchange $\Omega_{a}$ and $\Omega_{b}$ so that now
$r=1/20$. The $\sigma$ values are now closer to those from EVFCS,
$\sigma_{a}=2.06$ and $\sigma_{b}=2.06\times 10^{-3}$, but due to the
proportions of the subdomains, we see a remarkable increase in $\sigma$
variance — leading to large $\phi_{\text{UB}}$. The contrast with DVFCSLF
demonstrates the importance of topology of the composite material with respect
to the environment in determining $\phi(\sigma)$.

Figure 15: Geometry Disparate Volume-Fraction Concentric Squares with Heavy
Film, denoted in short as DVFCSHF: solution $\psi^{\prime 0}$ and final
(adaptively refined) finite element mesh. Note we may construct from
$\Omega^{2\mathrm{D}}$ associated extruded domains in $\mathbb{R}^{3}$.

##  6 Error Estimation

In this section we establish QoI bounds for the solution to the initial value
problem (11) – (13) and errors associated with QoI approximations introduced
in Definition 4.2: $\tilde{u}_{\text{avg}}^{1}$,
$\tilde{u}_{\text{avg}}^{2\text{P}}$, and $\tilde{u}_{\Delta}^{2}$. We shall
consider primarily two domains: SART-1 and SART-2. Both are instances of our
Small-Angle Right Triangle introduced in Proposition 5.6 ; SART-1 corresponds
to a relatively small angle (equivalent to SART in Section 5.4); SART-2
corresponds to an even smaller angle ($W=1/16$). These domains serve in the
numerical examples presented throughout this section.

Also in this section it shall be convenient to introduce a “slow” time scale:

######  Definition 6.1 (Slow-Time Representation).

Let us define

$\displaystyle T\equiv B\gamma t,\ T_{\text{final}}\equiv B\gamma
t_{\text{final}}$ (89)

and associated field and quantities of interest,

$\displaystyle U(x,T)\equiv u(x,t=T/(B\gamma)),$ (90)

and

$\displaystyle U_{\text{avg}}(T)$ $\displaystyle\equiv
u_{\text{avg}}(t=T/(B\gamma)),$ (91) $\displaystyle
U_{\partial{\text{avg}}}(T)$ $\displaystyle\equiv
u_{\partial{\text{avg}}}(t=T/(B\gamma)),$ (92) $\displaystyle U_{\Delta}(T)$
$\displaystyle\equiv u_{\Delta}(t=T/(B\gamma)),$ (93)

respectively. The arguments $(x,T)$ (respectively, $T$) are an abbreviated
form for $(x,T;B)$ or even $(x,T;B;\sigma,\kappa)$ (respectively, $(T;B)$ or
even $(T;B;\sigma,\kappa)$).

In a similar fashion,

$\displaystyle\widetilde{U}_{\text{avg}}^{1}(T)$
$\displaystyle\equiv\tilde{u}_{\text{avg}}^{1}(t=T/(B\gamma)),$ (94)
$\displaystyle\widetilde{U}_{\text{avg}}^{2\text{P}}(T)$
$\displaystyle\equiv\tilde{u}_{\text{avg}}^{2\text{P}}(t=T/(B\gamma)),\text{
and}$ (95) $\displaystyle\widetilde{U}_{\Delta}^{2}(T)$
$\displaystyle\equiv\tilde{u}_{\Delta}^{2}(t=T/(B\gamma))$ (96)

are the slow-time representations of our lumped approximations of Definition
4.2. Again, $T$ is an abbreviated argument list.

Note that $\widetilde{U}_{\text{avg}}^{1}(T)=\exp(-T)$, which demonstrates the
convenience of the slow time representation — in which time is effectively
normalized with respect to the small-$B$ time constant.

Finally, we make a remark on nomenclature: $T=B\gamma
t={\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{h|\partial\Omega|t/(\rho
c|\Omega|)}}$ is our nondimensional slow time; in contrast,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{T}}$ is
the dimensional temperature, with nondimensional counterpart
$u\equiv{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{(T-T_{\infty})/(T_{\text{i}}-T_{\infty})}}$.
We shall not invoke dimensional variables in the remainder of this section.

###  6.1 First-Order Lumped Approximation: Domain-Average

Before we proceed to the main domain-average analysis, we first consider the
positivity of the relative domain-boundary difference $U_{\Delta}(T)$ for all
$T>0$:

######  Lemma 6.2 (Strictly Positive $U_{\Delta}$).

The relative difference between the domain-mean and the boundary-mean
temperature satisfies

$\displaystyle U_{\Delta}(T;B)>0,\ \forall\,T,B>0.$ (97)

For $T=0$ or $B=0$, we have equality of the domain and boundary averages such
that

$\displaystyle U_{\Delta}(0;B)$ $\displaystyle=0,$ (98) $\displaystyle
U_{\Delta}(T;0)$ $\displaystyle=0.$ (99)

We provide a sketch of the proof for Lemma 6.2 in Appendix C.14. Following
this result, we have:

######  Proposition 6.3 (Global Domain-Average Lower Bound).

The first-order lumped approximation $\widetilde{U}_{\text{avg}}^{1}$ is a
lower bound for the domain-average $U_{\text{avg}}(T)$ for all
$T,B\in\mathbb{R}_{0+}$. More specifically, for $B>0$ and $T>0$, we have a
strict inequality

$\displaystyle U_{\text{avg}}(T)>\widetilde{U}_{\text{avg}}^{1}(T).$ (100)

We note that this lower bound result has a simple physical interpretation: the
(first-order) lumped approximation neglects the conduction resistance,
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\mathscr{R}_{\text{avg}}}}$
of (56) in Section 4, and hence the associated heat transfer rate is greater
than the actual heat transfer rate, and, as a result, the lumped temperature
approaches equilibrium faster than the actual temperature, and thus — given
our nondimensionalization — the nondimensional lumped temperature will always
be less than the actual temperature. We provide a sketch of the proof for
Proposition 6.3 in Appendix C.15.

Now We define our errors associated with the first-order lumped approximation
$\widetilde{U}_{\text{avg}}^{1}$:

######  Definition 6.4 (First-Order Lumped Domain-Average: Absolute Error).

We define

$\displaystyle
E_{\text{avg}}^{1}({T_{\text{final}}})\equiv\max_{T\in[0,T_{\text{final}}]}U_{\text{avg}}(T)-\widetilde{U}_{\text{avg}}^{1}(T),$
(101)

As always, $T$ is an abbreviated argument list.

We can then prove, in similar fashion to [4], but for general geometry,

######  Proposition 6.5 (Asymptotic Error Estimates for
$\widetilde{U}_{\text{avg}}^{1}$).

The absolute error, for any final time, satisfies

$\displaystyle E_{\text{avg}}^{1}({T_{\text{final}}})\leq\frac{\phi
B}{\gamma\,e}+\pazocal{O}(B^{2})\text{ as }B\rightarrow 0\,.$ (102)

We note the leading error term

$\displaystyle E_{\text{avg}}^{{1}\,\text{asymp}}\equiv\underbrace{\frac{\phi
B}{\gamma}}_{\text{Bi}^{\prime}_{\text{dunk}}}e^{-1}\,$ (103)

is proportional to the corrected intrinsic Biot number,
$\text{Bi}^{\prime}_{\text{dunk}}$. We also note that

$\displaystyle T_{\text{max}}=1+\dfrac{1}{2}\underbrace{\frac{\phi
B}{\gamma}}_{\text{Bi}^{\prime}_{\text{dunk}}}+\,\pazocal{O}(B^{2})\text{ as
}B\rightarrow 0$ (104)

is the (slow) time $T$ at which
$U_{\text{avg}}(T)-\widetilde{U}_{\text{avg}}^{1}(T)$ attains the maximum. For
both expansions about $B=0$, the high-order terms encapsulated in
$\pazocal{O}(B^{2})$ are independent of ${T_{\text{final}}}$.

We note that in practice slow times $T$ of order unity are of greatest
interest. We provide in Appendix C.16 a sketch of the full proof of
Proposition 6.5. We also note that the relative error associated with
$\widetilde{U}_{\text{avg}}^{1}$ based on Proposition 6.5 can be defined (as
well as its associated error bound).444 We can define
$E_{\text{avg}}^{{1}\,\text{rel}}(T)\equiv\left(U_{\text{avg}}(T)-\widetilde{U}_{\text{avg}}^{1}(T)\right)/\widetilde{U}_{\text{avg}}^{1}(T)$
as the relative error associated with the first-order lumped approximation
$\widetilde{U}_{\text{avg}}^{1}$. We then define the relative error bound
$E_{\text{avg}}^{{1}\,\text{rel
asymp}}(T)\equiv\text{Bi}^{\prime}_{\text{dunk}}T$ such that
$E_{\text{avg}}^{{1}\,\text{rel}}(T)\leq E_{\text{avg}}^{{1}\,\text{rel
asymp}}+\pazocal{O}(B^{2})\text{ as }B\rightarrow 0$.

Now we consider non-asymptotic error bounds for
$\widetilde{U}_{\text{avg}}^{1}$:

######  Proposition 6.6 (Strict Error Bounds for
$\widetilde{U}_{\text{avg}}^{1}$).

For all positive $B$, the absolute error satisfies

$\displaystyle
E_{\text{avg}}^{1}({T_{\text{final}}})\,\leq\,\frac{1}{2}\sqrt{\frac{\phi
B}{\gamma}}\,,\ \forall\,{T_{\text{final}}}\in[0,\infty).$ (105)

We note that for this non-asymptotic result, the error bound

$\displaystyle
E_{\text{avg}}^{{1}\,\text{UB}}\equiv\frac{1}{2}\biggl{(}\underbrace{\frac{\phi
B}{\gamma}}_{\text{Bi}^{\prime}_{\text{dunk}}}\biggr{)}^{1/2}\;$ (106)

also depends on the corrected intrinsic Biot number
$\text{Bi}^{\prime}_{\text{dunk}}$.

We note from Proposition 6.6 and Proposition 6.3 respectively that the domain
average can be bounded by $\widetilde{U}_{\text{avg}}^{1}$ and
$\widetilde{U}_{\text{avg}}^{1}+E_{\text{avg}}^{{1}\,\text{UB}}$. We provide a
sketch of the proof for Proposition 6.6 in Appendix C.17. Similarly to
Proposition 6.5, we can define a relative error bound for
$\widetilde{U}_{\text{avg}}^{1}$.555 We can define
$E_{\text{avg}}^{{1}\,\text{rel
UB}}(T)\equiv\frac{1}{2}\sqrt{\text{Bi}^{\prime}_{\text{dunk}}}\exp(T)$, where
$E_{\text{avg}}^{{1}\,\text{rel}}(T)\leq E_{\text{avg}}^{{1}\,\text{rel
UB}}(T)$.

The strict bound of Proposition 6.6 eliminates any uncertainty in the error
estimate. However, it is clear from Proposition 6.5 and the result for
$E_{\text{avg}}^{{1}\,\text{asymp}}$ that the strict bound
$E_{\text{avg}}^{{1}\,\text{UB}}$ is not sharp as $B\rightarrow 0$. On the
other hand, for large $B$, the strict bound $E_{\text{avg}}^{1\,\text{UB}}$
will often be irrelevant — greater than unity. Note given the maximum
principle, the error $E_{\text{avg}}^{1}$ must always be less than one.
Ironically, the strict bound $E_{\text{avg}}^{1\,\text{UB}}$ may be most
relevant for very small $B$, confirming without any doubt that the error is
less than some acceptable threshold.

We show in Table 5 (respectively, Table 6) the absolute error and associated
asymptotic error estimates and strict error bounds for the first-order lumped
approximation of $U_{\text{avg}}$ for domain SART-1 (respectively, SART-2).
Note when $E_{\text{avg}}^{{1}\,\text{UB}}>1$, no useful information is
obtained and we hence indicate — in the tables. Note for SART-1, $\phi=9.136$,
and for SART-2, $\phi=161.1$. The results of these tables are obtained by
finite difference in time (BDF1 for the first time-step and BDF2 for
subsequent time-steps) and finite element in space ($\mathbb{P}_{2}$)
discretizations of the Robin heat equation. We integrate to
$T_{\text{final}}=2$ and present $E^{1}_{\text{avg}}(2)$. For any given
result, we consider a sequence of uniform refinements in space and time until
an error estimator (the $H^{1}$ norm of the discretization error in $U$ at
$T=T_{\text{final}}=2$) is suitably small compared to the predicted error
$E_{\text{avg}}^{1}$.

We make the following observations:

  1. 1.

We confirm, in both Table 5 and Table 6, that $E_{\text{avg}}^{1}$ indeed
tends to zero linearly in B: the first-order approximation is indeed first-
order.

  2. 2.

We confirm — by comparison of Table 5 and Table 6 — that, for a fixed $B$, or
$\text{Bi}_{\text{dunk}}$, the error increases substantially with increased
$\phi$.

  3. 3.

We confirm that $E_{\text{avg}}^{{1}\,\text{UB}}$ is indeed an upper bound,
albeit — and as expected — not sharp as $B\rightarrow 0$.

  4. 4.

We observe that
$E_{\text{avg}}^{{1}\,\text{asymp}}/E_{\text{avg}}^{1}\rightarrow 1$ and from
above as $B\rightarrow 0$: the asymptotic error estimate is a sharp upper
bound for small $B$. It is crucial to note that, absent the correct value of
$\phi$, $E_{\text{avg}}^{{1}\,\text{asymp}}/E_{\text{avg}}^{1}$ would tend to
a constant less than unity — underestimation of the error — or a constant
greater than unity — overestimation of the error — as $B$ tends to zero.

We first note that the numerical results are illustrations of the provided
mathematical proofs with regard to these observations. We also note that if
$E_{\text{avg}}^{{1}\,\text{asymp}}$ is constructed in terms of an upper bound
for $\phi$, for example as in Section 5.1 Proposition 5.4, then
$E_{\text{avg}}^{{1}\,\text{asymp}}$ will still provide an upper bound for
$E_{\text{avg}}^{1}$ as $B\rightarrow 0$, however
$E_{\text{avg}}^{{1}\,\text{asymp}}/E_{\text{avg}}^{1}$ will approach a
constant greater than unity. We have not provide a rigorous proof for this
statement, however.

Table 5: Absolute error and associated asymptotic error estimates and strict error bounds for the first-order lumped approximation of $U_{\text{avg}}(T)$, $T\in[0,2]$, for domain SART-1. $B$ | $\text{Bi}_{\text{dunk}}\equiv B/\gamma$ | $\text{Bi}^{\prime}_{\text{dunk}}\equiv\phi B/\gamma$ | $E^{1}_{\text{avg}}(2)$ | $E_{\text{avg}}^{\text{1 asymp}}$ | $E_{\text{avg}}^{\text{1 asymp}}/E^{1}_{\text{avg}}(2)$ | $E_{\text{avg}}^{\text{1 UB}}$
---|---|---|---|---|---|---
$1\;\\!\cdot\;\\!10^{-3}$ | $5.48\;\\!\cdot\;\\!10^{-5}$ | $5.01\;\\!\cdot\;\\!10^{-4}$ | $1.84\;\\!\cdot\;\\!10^{-4}$ | $1.84\;\\!\cdot\;\\!10^{-4}$ | $1.00$ | $1.12\;\\!\cdot\;\\!10^{-2}$
$2\;\\!\cdot\;\\!10^{-3}$ | $1.10\;\\!\cdot\;\\!10^{-4}$ | $1.00\;\\!\cdot\;\\!10^{-3}$ | $3.67\;\\!\cdot\;\\!10^{-4}$ | $3.68\;\\!\cdot\;\\!10^{-4}$ | $1.00$ | $1.58\;\\!\cdot\;\\!10^{-2}$
$5\;\\!\cdot\;\\!10^{-3}$ | $2.74\;\\!\cdot\;\\!10^{-4}$ | $2.50\;\\!\cdot\;\\!10^{-3}$ | $9.10\;\\!\cdot\;\\!10^{-4}$ | $9.21\;\\!\cdot\;\\!10^{-4}$ | $1.01$ | $2.50\;\\!\cdot\;\\!10^{-2}$
$1\;\\!\cdot\;\\!10^{-2}$ | $5.48\;\\!\cdot\;\\!10^{-4}$ | $5.01\;\\!\cdot\;\\!10^{-3}$ | $1.80\;\\!\cdot\;\\!10^{-3}$ | $1.84\;\\!\cdot\;\\!10^{-3}$ | $1.03$ | $3.54\;\\!\cdot\;\\!10^{-2}$
$2\;\\!\cdot\;\\!10^{-2}$ | $1.10\;\\!\cdot\;\\!10^{-3}$ | $1.00\;\\!\cdot\;\\!10^{-2}$ | $3.51\;\\!\cdot\;\\!10^{-3}$ | $3.68\;\\!\cdot\;\\!10^{-3}$ | $1.05$ | $5.00\;\\!\cdot\;\\!10^{-2}$
$5\;\\!\cdot\;\\!10^{-2}$ | $2.74\;\\!\cdot\;\\!10^{-3}$ | $2.50\;\\!\cdot\;\\!10^{-2}$ | $8.17\;\\!\cdot\;\\!10^{-3}$ | $9.21\;\\!\cdot\;\\!10^{-3}$ | $1.13$ | $7.91\;\\!\cdot\;\\!10^{-2}$
$1\;\\!\cdot\;\\!10^{-1}$ | $5.48\;\\!\cdot\;\\!10^{-3}$ | $5.01\;\\!\cdot\;\\!10^{-2}$ | $1.47\;\\!\cdot\;\\!10^{-2}$ | $1.84\;\\!\cdot\;\\!10^{-2}$ | $1.26$ | $1.12\;\\!\cdot\;\\!10^{-1}$
$2\;\\!\cdot\;\\!10^{-1}$ | $1.10\;\\!\cdot\;\\!10^{-2}$ | $1.00\;\\!\cdot\;\\!10^{-1}$ | $2.44\;\\!\cdot\;\\!10^{-2}$ | $3.68\;\\!\cdot\;\\!10^{-2}$ | $1.51$ | $1.58\;\\!\cdot\;\\!10^{-1}$
$5\;\\!\cdot\;\\!10^{-1}$ | $2.74\;\\!\cdot\;\\!10^{-2}$ | $2.50\;\\!\cdot\;\\!10^{-1}$ | $4.10\;\\!\cdot\;\\!10^{-2}$ | $9.21\;\\!\cdot\;\\!10^{-2}$ | $2.25$ | $2.50\;\\!\cdot\;\\!10^{-1}$
$1\;\\!\cdot\;\\!10^{0}$ | $5.48\;\\!\cdot\;\\!10^{-2}$ | $5.01\;\\!\cdot\;\\!10^{-1}$ | $5.55\;\\!\cdot\;\\!10^{-2}$ | $1.84\;\\!\cdot\;\\!10^{-1}$ | $3.32$ | $3.54\;\\!\cdot\;\\!10^{-1}$
Table 6: Absolute error and associated asymptotic error estimates and strict error bounds for the first-order lumped approximation of $U_{\text{avg}}(T)$, $T\in[0,2]$, for domain SART-2. $B$ | $\text{Bi}_{\text{dunk}}\equiv B/\gamma$ | $\text{Bi}^{\prime}_{\text{dunk}}\equiv\phi B/\gamma$ | $E^{1}_{\text{avg}}(2)$ | $E_{\text{avg}}^{\text{1 asymp}}$ | $E_{\text{avg}}^{\text{1 asymp}}/E^{1}_{\text{avg}}(2)$ | $E_{\text{avg}}^{\text{1 UB}}$
---|---|---|---|---|---|---
$1\;\\!\cdot\;\\!10^{-3}$ | $1.51\;\\!\cdot\;\\!10^{-5}$ | $\quad 2.44\;\\!\cdot\;\\!10^{-3}$ | $8.89\;\\!\cdot\;\\!10^{-4}$ | $8.97\;\\!\cdot\;\\!10^{-4}$ | $1.01$ | $2.47\;\\!\cdot\;\\!10^{-2}$
$2\;\\!\cdot\;\\!10^{-3}$ | $3.03\;\\!\cdot\;\\!10^{-5}$ | $\quad 4.88\;\\!\cdot\;\\!10^{-3}$ | $1.76\;\\!\cdot\;\\!10^{-3}$ | $1.79\;\\!\cdot\;\\!10^{-3}$ | $1.02$ | $3.49\;\\!\cdot\;\\!10^{-2}$
$5\;\\!\cdot\;\\!10^{-3}$ | $7.57\;\\!\cdot\;\\!10^{-5}$ | $\quad 1.22\;\\!\cdot\;\\!10^{-2}$ | $4.27\;\\!\cdot\;\\!10^{-3}$ | $4.49\;\\!\cdot\;\\!10^{-3}$ | $1.05$ | $5.52\;\\!\cdot\;\\!10^{-2}$
$1\;\\!\cdot\;\\!10^{-2}$ | $1.51\;\\!\cdot\;\\!10^{-4}$ | $\quad 2.44\;\\!\cdot\;\\!10^{-2}$ | $8.14\;\\!\cdot\;\\!10^{-3}$ | $8.97\;\\!\cdot\;\\!10^{-3}$ | $1.10$ | $7.81\;\\!\cdot\;\\!10^{-2}$
$2\;\\!\cdot\;\\!10^{-2}$ | $3.03\;\\!\cdot\;\\!10^{-4}$ | $\quad 4.88\;\\!\cdot\;\\!10^{-2}$ | $1.49\;\\!\cdot\;\\!10^{-2}$ | $1.79\;\\!\cdot\;\\!10^{-2}$ | $1.21$ | $1.10\;\\!\cdot\;\\!10^{-1}$
$5\;\\!\cdot\;\\!10^{-2}$ | $7.57\;\\!\cdot\;\\!10^{-4}$ | $\quad 1.22\;\\!\cdot\;\\!10^{-1}$ | $2.93\;\\!\cdot\;\\!10^{-2}$ | $4.49\;\\!\cdot\;\\!10^{-2}$ | $1.53$ | $1.75\;\\!\cdot\;\\!10^{-1}$
$1\;\\!\cdot\;\\!10^{-1}$ | $1.51\;\\!\cdot\;\\!10^{-3}$ | $\quad 2.44\;\\!\cdot\;\\!10^{-1}$ | $4.31\;\\!\cdot\;\\!10^{-2}$ | $8.97\;\\!\cdot\;\\!10^{-2}$ | $2.08$ | $2.47\;\\!\cdot\;\\!10^{-1}$
$2\;\\!\cdot\;\\!10^{-1}$ | $3.03\;\\!\cdot\;\\!10^{-3}$ | $\quad 4.88\;\\!\cdot\;\\!10^{-1}$ | $5.59\;\\!\cdot\;\\!10^{-2}$ | $1.79\;\\!\cdot\;\\!10^{-1}$ | $3.21$ | $3.49\;\\!\cdot\;\\!10^{-1}$
$5\;\\!\cdot\;\\!10^{-1}$ | $7.57\;\\!\cdot\;\\!10^{-3}$ | $\quad 1.22\;\\!\cdot\;\\!10^{0}$ | $6.83\;\\!\cdot\;\\!10^{-2}$ | $4.49\;\\!\cdot\;\\!10^{-1}$ | $6.57$ | $5.52\;\\!\cdot\;\\!10^{-1}$
$1\;\\!\cdot\;\\!10^{0}$ | $1.51\;\\!\cdot\;\\!10^{-2}$ | $\quad 2.44\;\\!\cdot\;\\!10^{0}$ | $7.46\;\\!\cdot\;\\!10^{-2}$ | $8.97\;\\!\cdot\;\\!10^{-1}$ | $12.0$ | $7.81\;\\!\cdot\;\\!10^{-1}$

The approximation $\widetilde{U}_{\text{avg}}^{1}(T;B)$ is agnostic to the
details of material properties, but the approximation is not restricted to
problems characterized by uniform material properties. Note, however that
although the evaluation of the approximation is independent of $\sigma$ and
$\kappa$ — we require only $B$ and geometric information — the error estimate
must be defined relative to the (possibly) non-uniform problem data; the
latter appear in the equation for $\phi$.

###  6.2 Second-Order Lumped Approximation: Domain-Average

We now consider the second-order lumped approximation,
$\widetilde{U}_{\text{avg}}^{2\text{P}}$ of (95), which is based on a second-
order Padé-approximant (40) of $\lambda$. Note the Padé choice is preferred
over the Taylor series: in the former, unlike in the latter, $\lambda$ is
bounded as $B\rightarrow\infty$, with an attendant improvement in accuracy for
large $B$. Before proceeding, we introduce the third-order expansion for
$\lambda$:

######  Lemma 6.7 (First Eigenvalue Asymptotic Expansion Revisited).

The first eigenvalue admits expansion

$\displaystyle\lambda=\gamma B-\phi
B^{2}+(\chi-\gamma\,\Upsilon)B^{3}+\pazocal{O}(B^{4})\text{ as }B\rightarrow
0.$ (107)

We recall from (37) that
$\phi\equiv\int_{\Omega}\kappa\,\nabla\psi^{\prime\,0}\cdot\nabla\psi^{\prime\,0}$,
and from (35) that
$\Upsilon\equiv\int_{\Omega}\sigma\,\psi^{\prime\,0}\psi^{\prime\,0}$; we now
further introduce $\chi\in\mathbb{R}_{0+}$ given by

$\displaystyle\chi$
$\displaystyle\equiv\int_{\partial\Omega}\psi^{\prime\,0}\psi^{\prime\,0}.$
(108)

we recall the definition $\psi^{\prime\,0}\equiv d\psi_{1}/dB\text{ evaluated
at }B=0$.

We provide in Appendix C.18 a sketch of the proof of Proposition 6.7.

######  Proposition 6.8 (Scale-Invariance of $\gamma\chi$ and
$\gamma^{2}\Upsilon$).

We are given $\Omega_{1}\subset\mathbb{R}^{d}$ and a translation-rotation-
uniform dilation map of $\Omega_{1}$,
${\Omega}_{2}\equiv\\{\mathbb{M}(x)\,|\,x\in\Omega_{1}\,\\}$; here
$\mathbb{M}(x)=\mathbb{T}+\mathbb{Q}\,x+\alpha\,\mathbb{I}\,x$, where $\alpha$
is a positive real scalar, $\mathbb{T}\in\mathbb{R}^{d}$ is a translation
vector, $\mathbb{Q}\in\mathbb{R}^{d\times d}$ is a rotation matrix, and
$\mathbb{I}\in\mathbb{R}^{d\times d}$ is the identity. For given
$\sigma\equiv\sigma_{1},\kappa\equiv\kappa_{1}$, $\Omega\equiv\Omega_{1}$ we
evaluate $\chi_{1}\equiv\chi(\Omega_{1})$,
$\Upsilon_{1}\equiv\Upsilon(\Omega_{1})$, and
$\gamma_{1}\equiv|\Omega_{1}|/|\partial\Omega_{1}|$; for given
$\sigma\equiv\sigma_{2}\equiv\sigma_{1}\circ\mathbb{M}^{-1},\kappa\equiv\kappa_{2}\equiv\kappa_{1}\circ\mathbb{M}^{-1}$,
$\Omega\equiv\Omega_{2}$ we evaluate $\chi_{2}\equiv\chi(\Omega_{2})$,
$\Upsilon_{2}\equiv\Upsilon(\Omega_{2})$, and
$\gamma_{2}\equiv|\Omega_{2}|/|\partial\Omega_{2}|$. Then, for any $\alpha$,
$\mathbb{T}$, and $\mathbb{Q}$, $\gamma_{1}\,\chi_{1}=\gamma_{2}\,\chi_{2}$
and $\gamma_{1}^{2}\,\Upsilon_{1}=\gamma_{2}^{2}\,\Upsilon_{2}$.

We provide in Appendix C.19 a sketch of the proof of Proposition 6.8.

We now summarize the computation of $\phi$, $\Upsilon$, and $\chi$ for given
domain $\Omega$ and coefficients $\sigma$, $\kappa$: Step 1, we find
$\psi^{\prime\,0}_{h}(\approx\psi^{\prime\,0})$ by finite element
approximation of PDE (Appendix A.5); Step 2, we evaluate $\phi_{h}$,
$\Upsilon_{h}$, and $\chi_{h}$ – each of which is a quadratic functional of
$\psi^{\prime 0}$. Clearly, Step 1 dominates the computational cost — Step 2
is effectively free. It follows that the cost to compute $\Upsilon$ and $\chi$
— which shall be required to form our error estimator for
$\widetilde{U}_{\text{avg}}^{2\text{P}}$ — is effectively the same as the cost
to compute $\phi$ — already required to form
$\widetilde{U}_{\text{avg}}^{2\text{P}}$. It further follows that the marginal
computational cost associated with extension from first-order to second-order,
if we consider both prediction and error estimation, is negligible.

We now provide the error estimate associated with our second-order Padé
approximation:

######  Proposition 6.9 (Asymptotic Error Bound for
$\widetilde{U}_{\text{avg}}^{2\text{P}}$).

We define

$\displaystyle
E_{\text{avg}}^{2\text{P}}({T_{\text{final}}})\equiv\max_{T\in[0,T_{\text{final}}]}|U_{\text{avg}}(T)-\widetilde{U}_{\text{avg}}^{2\text{P}}(T)|,$
(109)

for $\widetilde{U}_{\text{avg}}^{2\text{P}}$ given by Definition 6.1 and the
original QoI, $\tilde{u}_{\text{avg}}^{2\text{P}}$, given by Definition 4.2.

Then the asymptotic error satisfies

$\displaystyle
E_{\text{avg}}^{2\text{P}}({T_{\text{final}}})=E_{\text{avg}}^{2\text{P\,asymp}}+\pazocal{O}(B^{3})\text{
as }B\rightarrow 0$ (110)

where

$\displaystyle
E_{\text{avg}}^{2\text{P\,asymp}}\equiv\left(\frac{|\gamma\chi-\gamma^{2}\Upsilon-\phi^{2}|}{e\gamma^{2}}+\Upsilon\right)B^{2};$
(111)

which may also be expressed with scale-invariant terms and
$\text{Bi}_{\text{dunk}}$ as

$\displaystyle
E_{\text{avg}}^{2\text{P\,asymp}}=\left(|\gamma\chi-\gamma^{2}\Upsilon-\phi^{2}|e^{-1}+\gamma^{2}\Upsilon\right)\text{Bi}_{\text{dunk}}^{2}.$
(112)

Recall that the second-order asymptotic bound
$E_{\text{avg}}^{2P\,\text{asymp}}$ requires no PDE solutions other than
$\psi^{\prime\,0}$; since $\phi$, $\chi$, and $\Upsilon$ are all functional
evaluations of $\psi^{\prime\,0}$. Note also that all scale-dependence in
(112) is reflected in $\text{Bi}^{\prime}_{\text{dunk}}$. We provide in
Appendix C.20 a sketch of the proof of Proposition 6.9.

Before proceeding to the numerical results for the second-order approximation,
we provide a table of normalized $\chi$ and $\Upsilon$ values for the
canonical domains and selected triangles in Table 7. Recall that $\gamma\chi$
and $\gamma^{2}\Upsilon$ are scale-invariant (Proposition 6.8). These values
are derived by evaluation of the appropriate quadratic functionals of
$\psi^{\prime\,0}$ on each of these domains (Appendix B for the case of
sphere, and Appendix C.11 for the selected triangles). Note for the selected
triangles, we must first solve for the constant $b_{00}$ to satisfy the
normalization constraint $M(\psi^{\prime\,0})=0$.

$\Omega$ | $\phi$ | $\gamma\,\chi$ | $\gamma^{2}\,\Upsilon$
---|---|---|---
$\Omega_{\text{ref,1}}$: Interval | 1/3 | 1/9 | 1/45
$\Omega_{\text{ref,2}}$: Disk | 1/2 | 1/4 | 1/12
$\Omega_{\text{ref,3}}$: Sphere | 3/5 | 9/25 | 27/175
$\Omega_{\text{right triangle}}(W=1)$ | 4/3 | $\frac{4}{5}(3+2\sqrt{2})$ | $\frac{4}{15}(3+2\sqrt{2})$
$\Omega_{\text{equilateral triangle}}$ | 1 | 9/5 | 3/5
Table 7: Values of normalized functional evaluations for canonical domains and
selected triangles.

Before proceeding to the error results for SART-1 and SART-2, we tabulate the
normalized functional evaluations for each geometry, necessary for error
estimation for $\widetilde{U}_{\text{avg}}^{2\text{P}}$, in Table 8. We
present the numerical results for the error $E_{\text{avg}}^{2\text{P}}(2)$
for the two domains in Tables 9 and 10, respectively. We observe similar
trends as for the first-order results; however we note that the bound for the
second-order approximation is looser than that of the first-order
approximation. This is due to uncertainty in the temporal behavior of the
higher eigenmodes: requiring an additional bounding term (244) which does not
decay in time. However, for low $B$, we see the ratio is less than 2.5. for
SART-1 and less than 2 for SART-2.

$\Omega$ | $\phi$ | $\gamma\,\chi$ | $\gamma^{2}\,\Upsilon$ | $e_{h}^{\phi}$
---|---|---|---|---
SART-1 | 9.14 | $4.65\;\\!\cdot\;\\!10^{2}$ | $1.55\;\\!\cdot\;\\!10^{2}$ | $2.07\;\\!\cdot\;\\!10^{-13}$
SART-2 | 161 | $1.21\;\\!\cdot\;\\!10^{5}$ | $4.02\;\\!\cdot\;\\!10^{4}$ | $1.57\;\\!\cdot\;\\!10^{-11}$
Table 8: Values of normalized $\chi$ and $\Upsilon$ for SART-1 and SART-2 Domains. $B$ | $\text{Bi}_{\text{dunk}}\equiv B/\gamma$ | $\text{Bi}^{\prime}_{\text{dunk}}\equiv\phi B/\gamma$ | $E_{\text{avg}}^{2\text{P}}(2)$ | $E_{\text{avg}}^{2\text{P\,asymp}}$ | $E_{\text{avg}}^{2\text{P\,asymp}}/E_{\text{avg}}^{2\text{P}}(2)$
---|---|---|---|---|---
$1\;\\!\cdot\;\\!10^{-3}$ | $5.48\;\\!\cdot\;\\!10^{-5}$ | $5.01\;\\!\cdot\;\\!10^{-4}$ | $3.24\;\\!\cdot\;\\!10^{-7}$ | $7.16\;\\!\cdot\;\\!10^{-7}$ | $2.21$
$2\;\\!\cdot\;\\!10^{-3}$ | $1.10\;\\!\cdot\;\\!10^{-4}$ | $1.00\;\\!\cdot\;\\!10^{-3}$ | $1.72\;\\!\cdot\;\\!10^{-6}$ | $2.86\;\\!\cdot\;\\!10^{-6}$ | $1.66$
$5\;\\!\cdot\;\\!10^{-3}$ | $2.74\;\\!\cdot\;\\!10^{-4}$ | $2.50\;\\!\cdot\;\\!10^{-3}$ | $1.19\;\\!\cdot\;\\!10^{-5}$ | $1.79\;\\!\cdot\;\\!10^{-5}$ | $1.50$
$1\;\\!\cdot\;\\!10^{-2}$ | $5.48\;\\!\cdot\;\\!10^{-4}$ | $5.01\;\\!\cdot\;\\!10^{-3}$ | $4.79\;\\!\cdot\;\\!10^{-5}$ | $7.16\;\\!\cdot\;\\!10^{-5}$ | $1.50$
$2\;\\!\cdot\;\\!10^{-2}$ | $1.10\;\\!\cdot\;\\!10^{-3}$ | $1.00\;\\!\cdot\;\\!10^{-2}$ | $1.87\;\\!\cdot\;\\!10^{-4}$ | $2.86\;\\!\cdot\;\\!10^{-4}$ | $1.53$
$5\;\\!\cdot\;\\!10^{-2}$ | $2.74\;\\!\cdot\;\\!10^{-3}$ | $2.50\;\\!\cdot\;\\!10^{-2}$ | $1.07\;\\!\cdot\;\\!10^{-3}$ | $1.79\;\\!\cdot\;\\!10^{-3}$ | $1.67$
$1\;\\!\cdot\;\\!10^{-1}$ | $5.48\;\\!\cdot\;\\!10^{-3}$ | $5.01\;\\!\cdot\;\\!10^{-2}$ | $3.73\;\\!\cdot\;\\!10^{-3}$ | $7.16\;\\!\cdot\;\\!10^{-3}$ | $1.92$
$2\;\\!\cdot\;\\!10^{-1}$ | $1.10\;\\!\cdot\;\\!10^{-2}$ | $1.00\;\\!\cdot\;\\!10^{-1}$ | $1.17\;\\!\cdot\;\\!10^{-2}$ | $2.86\;\\!\cdot\;\\!10^{-2}$ | $2.46$
$5\;\\!\cdot\;\\!10^{-1}$ | $2.74\;\\!\cdot\;\\!10^{-2}$ | $2.50\;\\!\cdot\;\\!10^{-1}$ | $4.23\;\\!\cdot\;\\!10^{-2}$ | $1.79\;\\!\cdot\;\\!10^{-1}$ | $4.24$
$1\;\\!\cdot\;\\!10^{0}$ | $5.48\;\\!\cdot\;\\!10^{-2}$ | $5.01\;\\!\cdot\;\\!10^{-1}$ | $9.35\;\\!\cdot\;\\!10^{-2}$ | $7.16\;\\!\cdot\;\\!10^{-1}$ | $7.66$
Table 9: Error in Approximation $\widetilde{U}_{\text{avg}}^{2}$ to $U_{\text{avg}}$ for $T_{\text{final}}=2$: SART-1 Domain. $B$ | $\text{Bi}_{\text{dunk}}\equiv B/\gamma$ | $\text{Bi}^{\prime}_{\text{dunk}}\equiv\phi B/\gamma$ | $E_{\text{avg}}^{2\text{P}}(2)$ | $E_{\text{avg}}^{2\text{P\,asymp}}$ | $E_{\text{avg}}^{2\text{P\,asymp}}/E_{\text{avg}}^{2\text{P}}(2)$
---|---|---|---|---|---
$1\;\\!\cdot\;\\!10^{-3}$ | $1.51\;\\!\cdot\;\\!10^{-5}$ | $\quad 2.44\;\\!\cdot\;\\!10^{-3}$ | $9.18\;\\!\cdot\;\\!10^{-6}$ | $1.38\;\\!\cdot\;\\!10^{-5}$ | $1.50$