A Bayesian Spatiotemporal Modeling Approach to the Inverse Heat Conduction Problem

Inverse heat conduction problems involve estimating one or more unknown parameters or boundary conditions within thermal systems or processes using indirect temperature, heat flux, or infrared measurements. In recent years, IHCP has found a multitude of applications, ranging from the design of rocket nozzles for spaceflights [1], to the characterization and testing of materials under extreme thermal conditions [2], the evaluation of thermal efficiency in industrial processes such as reheating furnaces [3], and its implementation in nuclear power plants [4]. Moreover, IHCP has ventured into the realm of medical science, with its recent utilization in breast cancer detection and tumor localization [5].

The earliest works on solving IHCP date back to Stolz [6] and Burggraf [7]. In his pioneering work, Stolz treated the inverse problem as though it were of the forward type by numerically integrating it for unknown surface conditions followed by a graphical method for inverting the integrals as described in Ref. [8]. Although not specifically mentioned, Stolz realized the challenges emerging from the ill-posedness of IHCP, particularly the sensitivity of experimental error on the inverse parameter estimation and provided empirical methods to remedy them. However, it was not until 1980 that Calderon provided the conditions for the uniqueness of the inverse problem [9]. Subsequent development of Tikhonov's regularization methods [10] and sequential function specification techniques [1] for uniquely solving inverse problems fueled the research on IHCP. Readers are referred to Refs. [1] and [11] for an in-depth review of the methods. It is crucial to note that a majority of the works that followed formulated the IHCP as a least squares problem over the whole domain followed by gradient optimization techniques [12–14]. However, these methods only provided a point estimate on the unknown quantity of interest, without any quantification of the uncertainty. Wang and Zabaras reformulated the inverse heat conduction as a Bayesian inference problem by considering the unknown quantity of interest as a random variable [15]. The resulting Bayesian methods provided a framework to incorporate the uncertainty in measurements into the inverse estimation of the unknown quantity of interests.

Traditionally, solving IHCP requires simulating the heat equation (also known as the forward problem). However, repeatedly solving the forward problem is time-intensive and demands considerable computational resources. An alternative involves the development of surrogate models aimed at approximating the solutions of the underlying governing equations. Meanwhile, it's noted that existing methods employing surrogate models for solving IHCP fall short of explicitly incorporating the spatiotemporal dynamics observed in temperature measurements obtained from sensors within the heat-conducting domain over time. Other examples of such spatiotemporal relationships are observed in electro-encephalogram source localization involving electrical signals recorded at different scalp locations [16] and fluid dynamics where spatiotemporally varying velocity fields are measured to estimate fluid properties [17]. Failure to account for spatiotemporal relationships may lead to inaccurate estimation of the quantity of interest. Existing methods often simplify solutions by ignoring the spatiotemporal structure [15,18] or treat time series as independent identically distributed observations [18,19]. Some scenarios only consider the time series property, neglecting spatial information, leading to insufficient modeling and suboptimal solutions.

This work proposes a Bayesian spatiotemporal approach for solving IHCPs. Our approach first considers a penalized spline, or a P-spline formulated within a spatiotemporal forward model that incorporates spatial, temporal, and interaction terms, as a surrogate model for the forward heat conduction problem. The proposed model captures the spatial and temporal variation of temperature together with the interaction effect based on a conditional autoregressive prior [20]. By accounting for the interaction effect between the spatial and temporal parameters, we ensure that the available spatiotemporal information in the data is completely utilized, leading to a more accurate forward model, thereby improving the performance of IHCP. Subsequently, we discuss a Gibbs sampling algorithm based on the Bayesian inference framework for the estimation and quantification of the uncertainty of the unknown temporal, spatial, and interaction parameters. We implement our approach for estimating the transient heat flux in a simulated example and a pool-boiling experiment using the temperature collected from the thermocouples in the domain. We summarize the main contributions of this work as follows.

1. In this work, we present a Bayesian spatiotemporal modeling and solution approach for IHCPs, utilizing penalized splines in a spatiotemporal forward model. The approach integrates spatial, temporal, and interaction terms, enhancing accuracy and solving the inverse problem through a Bayesian inference framework and Gibbs sampling for parameter estimation.

2. By applying this methodology to case studies, including a 1D heat simulation and a pool boiling experiment with thermocouple data, we demonstrate improved heat flux reconstruction. This illustrates the potential of the model beyond traditional methods, highlighting its applicability in solving IHCPs with complex spatiotemporal dynamics.

The remainder of this paper is structured as follows: Section 2 provides a brief review of the literature on IHCP. In Sec. 3, we present a brief description of the IHCP of interest along with the proposed spatiotemporal modeling approach. The case study and discussion are provided in Sec. 4, and the paper is concluded in Sec. 5.

Heat conduction problems can be divided into two main categories: direct and inverse problems. Direct heat conduction problems involve estimating the temperature profile within a conducting solid given known boundary and initial conditions, heat generation rates, and thermophysical properties [21]. In contrast, IHCPs aim to determine surface temperatures, heat source rates, heat flux, and other thermophysical properties using measured temperatures within the conducting medium [1]. Direct heat conduction problems are generally well-posed, whereas IHCPs are mathematically ill-posed [22]. This is due to the fact that unavoidable random errors (noise) in the measured data lead to errors in the estimated unknowns that increase by several orders of magnitude [1]. While solutions to the forward heat conduction problems are nontrivial, a significant amount of research has been directed into exploring methods that can provide unique and stable solutions to this problem more efficiently [23]. Therefore, standard solution methods for solving different variations and problem focus of the forward heat conduction problem have been proposed in the literature [15]. The inverse heat conduction problem on the other hand is an active area of research. This is because of the complexity of the IHCPs arising from ill-posedness from the embedded noise in the measured data. The categorization of methods for solving IHCPs can be approached from two different perspectives: the solution framework [24] and the modeling framework [21].

From the modeling framework perspective, two main solution approaches have been developed: the mesh approach and the meshless approach. Extensive research has been conducted on both solution approaches and has been presented in the literature. For example, over the last few decades, a variety of mesh-based schemes have been explored to solve IHCPs. These include space marching difference schemes [25], slowly divergent space marching schemes [26], and the mollification method [27]. Other approaches that have been investigated include the control volume algorithm combined with a digital filter method [28], the difference approximation method [29], and the use of the finite element method (FEM) [30], finite difference methods (FDM) and Fourier transform techniques [31]. Meshless approaches have gained significant attention in the past two decades. For example, the method of fundamental solution (MFS) has been explored [32], along with the Legendre polynomials method [33], the meshless local Petrov–Galerkin (MLPG) method [34], and the use of radial basis functions (RBFs) and the thin-plate spline (TPS) radial basis functions method [35]. Interested readers may consult [21] for a more extensive review of these methods.

From the perspective of solving the IHCPs, the methods can be categorized as deterministic or Bayesian. The deterministic methods encompass a broad category of solution methodology, including the regularization methods [10,35,36,37,38], computational methods [25,26,30,31], optimization methods [39,40], and deep learning methods [41,42]. On the other hand, some of the studies that have adopted the Bayesian solution framework for the IHCP include Refs. [15,43,24,44,45]. The Bayesian framework was first demonstrated as a viable approach for solving inverse problems that arise in geophysical applications in the late 1960s and early 1970s [46,47]. The Bayesian approach to inverse problems has been further developed and improved theoretically and in connection with applications in heat transfer and medical imaging, in the following decades, as summarized in Ref. [48]. For an extensive review of Bayesian methodology with different cases, see Ref. [49]. Comparing regularization methods with Bayesian inversion results can be challenging due to the comprehensive nature of Bayesian modeling. Bayesian approaches facilitate the incorporation of additional information regarding the problem of interest, such as the inclusion of domain knowledge through prior modeling, a feature not available in regularization methods. Additionally, while regularization methods aim to provide a single stable numerical solution, Bayesian inversion seeks to furnish point estimates along with the posterior uncertainty [24].

Broadly speaking, the Bayesian framework provides a straightforward way to incorporate data from different sources into a single formulation. This advantageous property of the Bayesian solution framework makes it a convenient option for spatiotemporal modeling of the IHCP. By incorporating information from both the spatial and temporal domains into the data, the framework enables improved learning of unknown parameters. Spatiotemporal models have proved practical for parameter estimation purposes involving space–time-varying data in several engineering and scientific domains. A diverse array of inverse problems exists whose observed measurement is concurrently sampled spatially and temporally, and in which the data-generating process emphasizes a connection between close spatial and temporal instances. These inverse problems are called spatiotemporal inverse problems, for instance, the inverse heat conduction problem. In a spatiotemporal inverse problem, observed multivariate time series are used to infer parameters of physical or phenological interest from predefined data models [50]. Existing methods for solving these inverse problems often overlook the spatiotemporal information contained in the data (static model) or rather define a data model that summarizes the data averaged over time [51]. In either case, the data information involving the spatiotemporal interactions is not fully utilized for parametric inference. This leads to inadequate modeling of these problems, the effect of which results in the suboptimal performance of the model during evaluation. Recently, integration of the Bayesian solution framework with spatiotemporal modeling has been carried out in different studies such as Refs. [20] and [51]. However, these approaches did not model the spatiotemporal interaction terms in estimating the unknown quantity of interest. This limits the performance of the model, as the information encoded in the interaction parameter is completely left unused, as seen in Ref. [11]. Surprisingly, explicit modeling of the interaction effect when defining model parameters has been reported in the literature to have led to improved model performance [52,53]. Therefore, we adopt the spatiotemporal modeling approach for this problem type that includes the spatiotemporal interaction effect in the model terms. Moreover, we use the Bayesian solution methodology to achieve a robust solution to the problem to estimate the posterior distributions of the parameters of interest.

In a heat conduction problem, thermal energy is supplied to a conducting material from an external source at some known region in domain Ω. The resulting heat flows through the conducting material and exits through the boundary ∂Ω. Thermocouples are attached to the conducting solid at specific locations in Ω to measure the temperatures over time. Spatially distributed thermocouples recording temperature over time generate a spatiotemporal structure. Figure 1 shows the schematic diagram of a conducting solid with different boundaries surrounding the domain Ω. As with any measurement system, the temperatures recorded at the thermocouples are noisy. The initial and boundary conditions are provided for some boundary regions ∂Ω_D and ∂Ω_N in the form of temperature measurements or heat flux along the boundary. In the following, we provide details of the forward and inverse heat conduction problem, followed by the proposed spatiotemporal model for solving this problem.

Here, κ is the thermal conductivity, $c=ρCp$ is the heat capacity, r represents the position vector, ρ is the density of conducting solid and C_p is the heat capacity at constant pressure. The initial condition Eq. (4) and boundary conditions Eqs. (2) and (3) include the initial temperature T₀(r) on the entire domain Ω, the boundary temperature T_D(t,r) at ∂Ω_D and the boundary heat flux $qh$(t,r) on the subdomain ∂Ω_D and ∂Ω_N, respectively (see Fig. 1). The heat source or sink term is denoted by q(t,r). In solving the forward problem, we denote the computed temperatures at each time instance t, t = 1, …, N, and sensor location $rj$⁠, j = 1, …, M as T (t, $rj$⁠; $qh$⁠). The inverse heat conduction problem is detailed in Sec. 3.2.

The inverse heat conduction problem aims to estimate unknown parameters, such as boundary heat flux q_h(t,r), or the thermal conductivity κ, given known conditions elsewhere within the conducting solid Ω and temperature measurements from the sensors. With the solution of the forward problem available from computational simulation of governing equations [31] or surrogate models [40] and temperature measurements from the thermocouples, the unknown quantity of interest is estimated. In this work, the unknown quantity of interest is the boundary heat flux. Therefore, the IHCP concerned in this work aims to estimate the unknown heat flux in some boundary within the domain by leveraging the forward heat conduction equations [Eqs. (1)–(4)] and noisy temperature measurements. Traditional approaches to IHCP require repeatedly solving the forward problem by simulating the heat equation as described in Eqs. (1)–(4). However, this is time-intensive and requires significant computational resources. Alternatively, surrogate models are built to approximate the solution of the governing equations. As mentioned earlier, existing surrogate models do not explicitly account for the spatiotemporal relationships that exist in the temperature recorded at one or more sensors inside the heat conducting domain over time. Our proposed solution approach overcomes this limitation by explicitly modeling the spatiotemporal effects to accurately capture the nonlinear variations of temperature in space and time, e.g., during critical heat flux (CHF), in estimating the heat flux.

where $Ytj$ represents the temperature recorded at time instance t and sensor location $rj, ωtj$ is the measurement noise, $T (t, rj; qh)$ is the true temperature at location $rj$ and time t for some heat flux $qh$⁠. In the IHCP, the solutions to the forward problem, denoted as $T (t, rj; qh)$ at all time instances and locations, allows for the estimation of $qh$⁠, such that the temperature predictions from the estimated heat flux match the measured values.

where ω represents the flattened residual noise vector defined as ω = [$ωtj$], t = 1, 2, …, N and j = 1, 2, …,$M.$ Given the structure of the forward model in Eq. (6), the inverse heat conduction problem transforms into a parameter estimation problem where the unknown quantities are θ, ϕ, and ψ. However, given the ill-posed nature of the inverse form of Eq. (6), a direct inversion for the estimation of θ, ϕ, and ψ, is not feasible. Bayesian inferences handle ill-posedness by (a) imposing a prior on the space of possible solutions of unknown parameters, thereby penalizing less likely values, (b) accounting for measurement errors into the likelihood function, and (c) propagating the uncertainty in measurements to the posterior distribution. The statistical distribution of the noise, often assumed to be Gaussian, is incorporated into the model. When the noise parameter is considered unknown, augmented forms are defined for both temporal and spatial parameters, and the augmented vectors are inferred from the measurement vector Y. This study focuses on the formulation of F (θ, ϕ, ψ)) in Eq. (6), which is modeled using a spatio-temporal framework as elaborated in Sec. 3.3. The Bayesian inference approach is utilized to estimate the marginal posterior distribution of the model parameters, with a Gibbs sampling algorithm designed for this purpose.

where t represents the time vector, and r is the vector of thermocouple locations within the domain at which the temperature is measured. T_I is the solution to the forward problem, with the zero-heat flux condition on ∂Ω_h (the limit at which the heat flux is to be computed) and the known initial conditions in ∂Ω_D and the boundary conditions in ∂Ω_N of the conducting domain. $f(t, θ)$ and $g(r, ϕ)$ capture the temporal and spatial variations in temperature. These functions are expressed in terms of triangular or hat basis matrices (see, e.g., Fig. 2), and their corresponding temporal and spatial parameters, which will be explained later. Specifically, $f(t, θ)$captures the temporal evolution of the temperature profile. Similarly, $g(r, ϕ)$ learns the spatial temperature distribution across space. Finally, the spatiotemporal interaction function is denoted by $h(t, r, ψ)$⁠, and $ω$ represents the error of the residual measurements.

The temporal function $f(t, θ)$ is defined as $f(t, θ)=(Bt⊗1M )θ,$ where $Bt$ is a basis matrix of dimension $N×m$ constructed using triangular or hat basis functions [54,55] with $m$ nodes evaluated at each time instant. Intuitively, $m$ controls the resolution of the temporal discretization, which, in turn, directly influences the precision of the derived solutions in time. Figure 2(a) shows the basis function used in constructing this matrix. $1M$ is a column vector of ones of dimension $M×1$⁠, with M representing the number of thermocouples, and the ⊗ symbol represents the Kronecker product between $Bt$ and $1M.$Therefore, the dimension of the output vector from $f(t, θ)$ is therefore $n×1$⁠, where $n=NM$ is the total number of measurements recorded from all the thermocouples. The parameter $θ =(θ1, …, θm)′$ contains the coefficients of the temporal basis defined for all time instances.

where V = D′D and D is a difference matrix of order 1 and dimension (p − 1) × p, λ₂ is the regularization parameter for the spatial effect. This parameter is also selected through sensitivity analysis following the same approach described earlier for selecting λ₁.

Finally, the spatiotemporal interaction function h(t, r, ψ) is defined as h(t, r, ψ) = $((G2⊗1N) ⊙ (1M⊗G1))ψ$⁠, where G₁ and G₂ are N × N and M × M basis matrix generated using triangular basis function, as shown in Fig. 2 with N and M number of nodes, respectively. 1_N and 1_M are column vectors of ones of dimensions N × 1 and M × 1, respectively. The symbol ⊙ is the element-wise matrix product. Given that the number of nodes in the basis is equal to the number of measurement steps N and thermocouples M for G₁ and G₂, respectively, the matrix obtained from the operation [(G₂$⊗ 1N$⁠) ⊙ (⁠$1M$ ⊗ G₁)] is thus an identity matrix of dimension n × n. This would not be the case if the number of nodes in the basis had been different from N and $M.$

where $λ3$ is a regularization parameter controlling the smoothness of the spatiotemporal effect. The selection of $λ3$ follows the same approach described for selecting $λ1$ and $λ2$⁠.

The residual error $ω$ is modeled as an additive white Gaussian noise with zero mean and variance $σ2$⁠. The variance parameter $σ2$ can be treated as an additional unknown to be inferred along with the other model parameters. However, to simplify the model, we estimate $σ2$ using sensitivity analysis as discussed in Sec. 3.4.1.

where $q(r, t)$ is the inferred heat flux, $Ŷ(r,t)$ are modeled probabilistic temperatures, and $κ$ is the known thermal conductivity. The calibrated spatiotemporal model enables uncertainty quantification in $Ŷ(r,t)$⁠, propagating to the estimated $qh.$

In the Sec. 3.3, we outlined the spatiotemporal model employed in this study, detailing the individual terms, their functions, and the parameters yet to be determined. The subsequent phase entails estimating these unknown parameters through the Bayesian inference framework.

where $p(θ, ϕ, ψ|Y)$is the joint posterior distribution of the temporal, spatial, and the interaction parameters that will be estimated, and $p(Y |θ, ϕ, ψ)$ denotes the likelihood. $p(θ), p(ϕ), and p(ψ)$ are the prior distributions of the temporal, spatial, and the interaction parameters, respectively.

where Y represents the vector of temperature measurements across time, from all sensors, with $Ytj$denoting the temperature measurement at a specific time instance t = 1 to N and thermocouple location j = 1 to $M. Ŷ$⁠, as defined in Eq. (7), is the model's predicted temperature for a given θ, ϕ, and ψ, with $Ŷtj$ corresponding to the predicted temperature at the same indices.

In the likelihood function, the noise variance $σ2$ encapsulates the variance of the residual errors between the model's predictions and the actual measurements. We consider two potential strategies for addressing σ²:

1. Treating σ² as an additional unknown parameter to be inferred alongside $θ, ϕ, and ψ.$ This approach allows the data to directly influence the determination of the noise variance, albeit at the cost of increased parameter space dimensionality.

2. Employing sensitivity analysis and leveraging prior knowledge to heuristically predetermine a reasonable value for $σ2.$ This strategy circumvents the expansion of the parameter space by an additional unknown.

In this study, we opt for the latter method to prevent the overparameterization of the model. Through sensitivity analyses concerning the levels of measurement noise, we determine a suitable value for σ² that adequately accounts for the anticipated residual errors. Although deriving σ² directly from the data could potentially enhance model accuracy, our preliminary analyses indicate that the model delivers satisfactory predictions across various noise levels. Consequently, a well-considered estimate of σ² suffices for our purposes, negating the need for its direct inference from the data. This fixed value of σ² also contributes to the efficiency of the Markov Chain Monte Carlo (MCMC) procedure by diminishing the quantity of unknown parameters. Future investigations will explore estimating σ² based on physical noise models. But for this initial study, using sensitivity analysis to select an appropriate noise variance a priori provides a reasonable model simplification that still captures residual errors.

After the posterior distribution of the unknown parameters has been formulated, the statistical analysis of the parameters can be carried out by exploring the posterior space. One of the most commonly used MCMC algorithms is the Gibbs sampling algorithm. According to Hoff [57], the Gibbs sampling algorithm is exceptionally adaptable for a multivariate distribution when the standard form of the conditional distribution of the individual random variable is known and can be sampled from. When the full conditional pdf is unknown, a proposal pdf is provided, and the drawn samples are either accepted or rejected based on an acceptable probability threshold. In situations where the full conditional pdfs are known, as it is for the current study where the conditional pdfs of $θ, ϕ, and ψ$ are known, these pdfs are used as the proposal distributions, and all drawn samples are accepted with a probability of 1. For the IHCP examined in this study, evaluating the likelihood function of Eq. (13) for each parameter set poses computational challenges, as it requires running numerical simulations to generate predicted temperatures. To address this, we developed a tailored Gibbs sampling algorithm leveraging the conditional distributions for our spatio-temporal model. In particular, we derive the full conditional posterior distributions for the parameter sets θ, ϕ, and ψ. Sampling directly from these conditionals enables efficient drawings from the joint posterior without needing to numerically evaluate the likelihood at every iteration. Our algorithm cycles through conditional sampling of each parameter in turn, exploiting the analytical conditional distributions to sidestep likelihood evaluations. This tailored Gibbs approach provides computational savings by avoiding explicit likelihood calculations during posterior exploration. The proposed Gibbs sampling algorithm is summarized below:

Initialize $θ(0), ϕ(0), ψ(0)$

For $i=0:Nmcmc$ (total sampled points)

Do:

• Sample $θj(i+1)∼ p (θj(i)θ−j(i), ϕ(i),ψ(i))$ for $j∈[1, m]$

• Sample $ϕk(i+1)∼ p (ϕk(i)ϕ−k(i), θ(i+1),ψ(i))$ for $k ∈[1, p]$

• Sample $ψl(i+1)∼ p (ψl(i)ψ−l(i), θ(i+1),ϕ(i+1))$ for $l∈[1, n]$

Return ${θ}i=1Nmcmc$⁠, ${ϕ}i=1Nmcmc$⁠, and ${ψ}i=1Nmcmc$

In this section, we showcase comprehensive applications of our proposed solution framework, which includes the spatiotemporal model and Gibbs sampling algorithm, by investigating two specific Inverse Heat Conduction Problems—a one-dimensional simulation study and a pool boiling experiment. We employ observed time series temperature data, collected via precision thermocouples at various points in the conducting region. The rich spatiotemporal characteristics of these data make it ideally suited for analysis using Spatio-Temporal Model (STM). Our primary objective is to reconstruct the posterior distribution of heat flux at a given boundary in the conductive material, utilizing the observed temperature data.

We assess the performance of the STM in comparison with three alternative approaches:

1. The Bayesian inverse model by Wang and Zabaras (W&Z) [15].

2. A spatiotemporal (ST) model lacking an interaction term.

3. A standard linear regression model described in Appendix  A.

By juxtaposing the predictions generated by the ST model with those from the STM, we highlight the pivotal role of the interaction term in the model. The significance of this interaction effect is systematically demonstrated across both case studies.

The properties of the conducting solid, used in the simulation, include thermal conductivity κ = 2 W/m K, density ρ = 10 kg/m³, and specific heat capacity C_p = 10 J/kg K. As detailed in Sec. 3.3, the parameters λ₁, λ₂, λ₃, and σ were selected using sensitivity analysis. Subsequently, we conducted this analysis and chose the values λ₁ = λ₂ = λ₃ = 1000, and σ = 0.01 for the current study. The outcomes highlight the superior efficacy of the spatiotemporal model compared to traditional methodologies under the outlined conditions. The reconstructed heat flux for each method is depicted in Fig. 3. The regression model exhibited the least accuracy in heat flux estimation, whereas our proposed model delivered optimal recovery across all methodologies. This assertion is corroborated in Table 1 within case study I, where the regression method registers the highest root-mean-squared-error (RMSE), and the STM method boasts the lowest RMSE. Concerning temperature prediction capabilities, each method's distinct forward model was employed to forecast the temperature based on their deduced heat flux. Spatial and temporal temperature predictions are showcased in Figs. 4 and 5, respectively. The corresponding RMSE values, calculated for each method, are enumerated in Table 1 under Case Study I. Yet again, the STM model outperforms in predicting temperatures both spatially and temporally, emphasizing its pre-eminence in heat flux recovery and temperature predictions over competing methods.

In this section, we applied the proposed solution framework, the developed STM, and the designed Gibb's sampling algorithm to solve a typical IHCP, the pool boiling experiment problem. This study uses experimental data from a pool boiling experiment for analysis. Pool boiling is a heat transfer mechanism that carries a phase transition from liquid to vapor [58].

In the pool boiling experiment, heat is supplied from a heating surface to a body of liquid (water) in contact with the surface. The setup of the pool boiling experiment is presented in Fig. 6 and a detailed description of the experimental method and conditions for the pool boiling experiment is provided in Appendix  B.

In this experiment, heat is supplied at the left boundary of the conducting medium, representing the source point. The heat is allowed to flow through a copper conductor to which four thermocouples have been attached before coming into contact with the pool at the right boundary. The time of study for this experiment is 556 s, and the sampling rate of each sensor Δt = 2s. The temporal discretization adopted for the heat flux q₀ is dt = 4s. Therefore, the values of $N, M, m, p$⁠, and n are 279, 4, 140, 5, and 1116, respectively. The measured temperatures from all the thermocouples are aggregated in vector $Y$⁠. We applied the proposed STM solution framework for heat flux reconstruction and temperature predictions and compare its performance with the conventional methods described earlier. The result is presented in Figs. 7 and 8 for the spatial and temporal temperature predictions, respectively. Using the RMSE values (shown in Table 2) of the competing methods as a performance index, the advantage of the proposed model is demonstrated.

In the study, we presented a new method for solving inverse heat conduction problems based on the spatiotemporal modeling approach and Bayesian inference framework. Modeling the spatiotemporal variations in temperature is critical to accurately estimating the unknown heat flux. We explicitly model the spatial, temporal, and interaction effect using a Gaussian random field prior. Finally, we proposed an MCMC sampling approach to draw samples from the posterior distribution of the unknown parameters. We benchmark the performance of the proposed methodology with other conventional methods using a 1D-simulation study and a pool boiling experimental case study. From the case studies, we demonstrate that the proposed method is superior in heat flux reconstruction and temperature predictions compared to other methods using standard metrics. We also highlight the importance of including the interaction terms in the model as demonstrated by the improved model performance. This result shows that the method has a potential to produce robust solutions to the inverse heat conduction problem. A future extension of this work will involve expanding the spatiotemporal model to handle three-dimensional transient heat conduction problems.

$Bs$ =

temporal basis matrix

$Bt$ =

spatial basis matrix

c =

heat capacity

C_p =

heat capacity at constant pressure FEM Finite Element Method

f =

temporal function

F =

direct system solver

g =

spatial function

h =

interaction function

I =

identity matrix

IHCP =

inverse heat conduction problems

m =

dimension of temporal parameter

M =

number of thermocouples

MCMC =

Markov Chain Monte Carlo

n =

total number of measurements

N =

number of measurement steps

p =

dimension of spatial parameter

P =

probability density function (PDF)

q =

deterministic inverse solution

q =

heat flux

q₀ =

known boundary heat flux

$qh$ =

unknown boundary heat flux

r =

thermocouple location

t =

time

T =

temperature

T₀ =

known initial temperature

T_D =

known boundary temperature

T_I =

initial and boundary temperature measurements

u =

random number

U =

covariance matrix of the interaction

V =

covariance matrix of the spatial MRF

w =

basis function of inverse solution

W =

covariance matrix of the temporal MRF

Y =

temperature measurement vector

∂Ω =

boundary region of the conducting solid

Δt =

thermocouple sampling interval

κ =

thermal conductivity

ρ =

density

Ω =

conducting region of the material

$δt$ =

time interval in discretization of the inverse solution

Given $x(t)=(x1, x2, …, xm)$ as the vector of the location of the thermocouples, and $t = (t0, t1, …, tn)$ as the discretized vector of time steps where the temperature has been sampled, a matrix of temperature measurements $T ∈ Rn×m$ is constructed containing sampled temperatures at each time-step. At every time instant $i ∈ {0, 1, …, n}$⁠, we solve a simple linear regression model with $x(t)$ as the feature vector and $T[i]$ as the response. The parameters (β0, β1) estimated at each time-step is then collected into the parameter matrix $B∈Rn×2$⁠. To estimate the temperatures at any location x in the domain across all time intervals, we use $T̂$ = B · x^⊺, where x^⊺ = [1, x] is the augmented location vector. We also estimate the heat flux as $q̂$ = −kβ1, where k is the thermal conductivity of the conducting region and β1 is the estimated slope coefficient.

The pool boiling test facility consists of a closed boiling chamber with a condensing loop for pressure control, a heating element with proper insulation and power supply, a multimodal sensing system, including temperature, pressure, and acoustic emission, and the associated data acquisition system. The boiling chamber consists of a polycarbonate ceiling and side walls, and a PEEK base. The ceiling of the boiling chamber is connected to a Graham condenser (Ace glass 5953-106) for degassing and a coiled copper tube condenser for pressure control, both of which are connected to a chiller (Thermo Scientific Polar ACCEL 500 Low/EA) using a bypass loop. Two compact screw-plug immersion heaters (McMaster-Carr 4668T54) are installed inside the boiling chamber for degassing and auxiliary heating. A variable-voltage transformer (McMaster-Carr 6994K17) is used to supply heat to the immersion heaters. The heating element includes a copper block with a boiling surface area of 10 × 10 mm², a custom-built PEEK insulation chamber and nine cartridge heaters (Omega Engineering HDC19102) inserted into the bottom of the copper block. The copper block is sealed using RTV (Momentive RTC 106) and epoxy (3M DP 110). The gap between the copper block and the PEEK chamber is filled with fiberglass and a ceramic block is used at the bottom of the PEAK chamber to support the copper block. The cartridge heaters are connected to a DC power supply (Magna-Power SL200-7.5) for heating power control. Four T-type thermocouples (Omega Engineering TJ36-CPSS-032 U-6) are evenly spaced in the copper block for surface temperature and heat flux measurements. A high-accuracy pressure transducer (Omega PX409-030A5V) is used to measure the vapor pressure inside the boiling chamber. Two T-type threaded pool thermocouples (McMaster-Carr 1245N16) are placed inside the boiling chamber to measure liquid and vapor temperatures, respectively. The signals from the thermocouples and pressure transducer are connected to an NI DAQ system (chassis: cDAQ-9178, temperature module: NI 9210, analogy input module: NI 9239), and recorded using LabVIEW. The experimental data used in this paper are obtained from steady-state pool boiling tests using de-ionized water. At the beginning of each experiment, water is degassed for half an hour with the Graham condenser open. After degassing, the Graham condenser is turned off and the coiled copper tube condenser is turned on for pressure control. At each heating power, the system reaches a steady-state in approximately 15 min and data sampling is performed during this time.