Reduced-Order Modeling of Pennes’ Bioheat Equation for Thermal Dose Analysis

Editor’s Note: The paper on which this article is based was originally presented at the 2023 IEEE International Symposium on Product Compliance Engineering (ISPCE), held in Dallas, TX in May 2023. It is reprinted here with the gracious permission of the IEEE. Copyright 2023, IEEE.

Introduction

Present users of wearable devices expect to be able to wear their devices for an entire day to make use of “always-on” functionalities such as step counting, heart-rate monitoring, and sleep quality tracking. This is especially true of wrist wearable devices such as smart watches and fitness trackers. In 2021, a record number of wearable devices were shipped for retail (over 530 million devices), representing a 20% increase in the number of devices shipped in 2020 [1]. As a result, more humans are maintaining longer duration contact with their powered electronic devices than ever before. Due to the necessity of direct contact between a user and a wearable device, some of the heat dissipated by the active elements of the device is transferred to the user’s skin. The objective of this study is to establish a methodology for building tractable “reduced order” models that can characterize and forecast the heat transfer between a device and a user and the potential for thermal injury. These reduced-order models can be implemented and solved efficiently in the context, for example, a control algorithm.

For the purposes of this study, heat that is dissipated from a powered device into the body is assumed to be conducted through four distinct tissue layers, as shown in Figure 1. The outermost layer is the epidermis, which forms a thin, protective barrier for the deeper layers of tissue. The epidermis is constantly shedding old cells but is maintained by new cell growth at the basal layer, which is the deepest portion of the epidermis. Beneath the basal layer is the dermis, which is a thicker layer of connective tissue that forms the core of the skin. The subcutaneous fat (also called the hypodermis) layer connects the dermis to the underlying inner tissue, composed of e.g., muscle and bone. Each of these tissues is distinct in their geometric, thermal, and physiological properties. In particular, the epidermis is not perfused by blood, while the other three tissues have vasculature. The perfusing blood can either remove or supply heat to these tissues, depending on their temperature. Skin burns are characterized by the depth of tissue that is damaged: first-degree burns damage but do not cause complete necrosis of the epidermis; second-degree burns cause complete necrosis of the epidermis but do not cause permanent damage to the dermis; and third-degree burns cause complete necrosis of the epidermis and significantly damage at least 75% of the dermis [2].

Quantitative thermal damage assessments are based on the temperatures experienced by the tissues and the duration of the thermal exposure. The thermal dose is often estimated in terms of Cumulative Equivalent Minutes at 43°C (CEM43°C) [3].

Any time-temperature history can be converted to an equivalent duration exposure at 43°C as follows:

    Eq. 1

where CEM43°C is the cumulative equivalent minutes at 43°C, T is the temperature of the tissue, t is physical time, and R is a piecewise-constant function of temperature given by:

    Eq. 2

Threshold values of CEM43°C that can lead to skin tissue damage are available from the literature, many of which have been derived from the work of Henriques and Moritz [4]. In the case of the skin tissues, the CEM43°C value that delineates injurious from non-injurious conditions ranges from approximately 300 to 600 minutes.

Governing Equations

The heat equation for conduction through a material containing a volumetric heat source, Q, is a partial differential equation (PDE) given by:

    Eq. 3

where ρ is the density of the material, c_p is the specific heat of the material, T is the temperature, t is physical time and k is the thermal conductivity of the material. The equation governing the thermal response of a tissue perfused by blood is instead modeled using Pennes’ Bioheat Equation as follows [5]:

    Eq. 4

where β is the “perfusion coefficient” (defined as the product of perfusion rate, blood density, and blood heat capacity), T_b is the temperature of the perfusing blood, and Q_met is the metabolic heat generation within the tissue. Note that (4) simplifies to (3) in the case that β = 0 and Q_met = Q. These equations can be solved by classical “full space” methods such as the finite difference or finite volume methods; however, these numerical methods become computationally expensive as the sizes of a simulation’s temporal and spatial domains grow. In particular, many modern consumer electronics contain a large number of different components and have complex internal geometries.

Proper Orthogonal Decomposition

Proper orthogonal decomposition (POD) is a numerical technique for projecting high-dimensional data into a lower-dimensional space that has been widely applied to heat, mass, and momentum transport problems. Given a dataset (e.g., temperature-time data from a thermal simulation), the POD technique finds a set of orthogonal basis vectors that will maximize the total variance of the data when projected onto these coordinates. In the parlance of POD, these orthogonal basis vectors are known as the spatial modes of the system.

The central assumption of the reduced-order model is that the temperature field calculated by the full space model can be approximated sufficiently well by a truncated sum of the first m POD spatial modes if they are appropriately modulated through time by temporal coefficients. This approximation can be written as:

    Eq. 5

where a_i(t) is the ith temporal coefficient, ϕ_i (x, y, z) is the ith spatial mode, and U₀(x, y, z) is a term that can be included to homogenize or center the temperature field. In particular, the inclusion of the U₀(x, y, z) term is useful for handling non-homogeneous Dirichlet (e.g., specified temperature) boundary conditions.

As suggested by (5), the temperature data from a thermal simulation may initially be stored in the form of an array with entries T (i, j, k, l), with the first index running over the simulation timesteps and the remaining indices running over the spatial coordinates. In order to use the POD technique, this array must first be reshaped into a two-dimensional matrix, T, in which each row contains the temperature at every point in the domain for a particular timestep. Once the data is organized in this way, the POD can be computed by several methods, including Singular Value Decomposition (SVD) or the Method of Snapshots [6]. The Method of Snapshots is the more efficient procedure when the number of timesteps (rows of T) is significantly smaller than the total number of spatial elements (columns of T), which is typically the case for simulation data. In the Method of Snapshots, the correlation matrix, χ, is first formed:

    Eq. 6

where N_t is the number of snapshots. Note that the size of the matrix TT^T is N_t × N_t and is therefore independent of the number of spatial elements in the simulation. An eigenvalue problem is then solved for the eigenvectors A and the matrix Ω with the corresponding eigenvalues on its diagonal:

    Eq. 7

The columns of A hold the values of each temporal coefficient at all snapshotted times. A matrix (denoted by Φ) containing the spatial modes of the system in its columns can be obtained from Φ = T^T A. Typically, each column of Φ is then normalized to unit magnitude, in which case the adjusted temporal coefficient values corresponding to the normalized spatial modes are given by the product of T and the normalized spatial mode matrix.

The number of POD modes retained in the approximation given by (5) influences the computational requirements of the reduced-order model and the fidelity of the approximation of the physics of the full space problem. Selecting too few POD modes may limit the ability of the reduced-order model to represent true physical phenomena, while retaining too many can introduce spurious dynamics or instability into the reduced-order model. A common heuristic is to retain modes such that a large fraction of the total variance of the entire mode spectrum is retained (this is sometimes referred to as the “energy” of the POD modes), i.e., choose m such that:

    Eq. 8

where λι is the ith eigenvalue of χ when ordered from largest to smallest and ε_tol is a user-defined tolerance.

Performing POD on a set of simulation data yields both spatial modes and temporal coefficients; however, these temporal coefficients only represent the dynamics in the original simulation data on which the decomposition was performed. Therefore, the temporal coefficients obtained during the decomposition are not useful in general. Instead, the temporal coefficients corresponding to arbitrary system dynamics can be calculated by projecting the full model onto a subset of the spatial modes, as detailed in the following section. The projection onto m spatial modes reduces the problem of solving the PDEs (3) and (4) to solving a system of m ordinary differential equations for the temporal coefficients.

Projection Methods

Once a reduced basis has been obtained from POD, the full physical model must be projected onto this low-dimensional space to form the ROM. There are two classes of methods that can be used to perform this projection: intrusive and non-intrusive methods. Intrusive methods require the modeler to know the differential and algebraic operators of the full space model and explicitly form their counterparts in the reduced space. In non-intrusive methods, the operators for the reduced order model are inferred directly from the same simulation data that was decomposed to find the POD modes.

Intrusive Projection of Pennes’ Bioheat Equation

The most common intrusive projection method used with POD-based ROMs is Galerkin’s method of weighted residuals. In this method, the approximation of the temperature field given by (5) is first substituted into the governing equations. This approximation introduces some error into the model. The method of weighted residuals attempts to minimize this error by requiring that the inner product of the model’s residual and a set of weight functions is equal to zero. When using the Galerkin method, these weight functions are taken to be the spatial modes of the system. Applying this methodology, with m spatial modes the original PDE (4) becomes a system of m ODEs:

    Eq. 9

where is the vector of temporal POD coefficients,  is the vector of the time derivatives of a, and the other matrices  and vector have elements:

where Ω represents the spatial domain of the simulation, Γ_D represents the boundary of Ω on which Dirichlet conditions are specified, and Γ_N represents the boundary of Ω on which Neumann conditions are specified (the boundary heat flux is denoted here by q). Note that the heat conduction term was integrated by parts to avoid the need for calculating second derivatives of numerical data. The initial condition for (9) is found by projecting the initial temperature condition T₀(x, y, z) ≡ T (t = 0, x, y, z) onto the spatial modes:

    Eq. 10

The matrices M and B and several terms in the vector f are not time-dependent and can therefore be pre-calculated prior to integrating the ODE system to speed up the solution process. The projection of the heat equation for non-perfuse tissue can also be recovered here by taking β equal to zero and replacing Q_met with an arbitrary (possibly time-varying) heat source.

Non-intrusive Projection of Pennes’ Bioheat Equation

Noting that Pennes’ Bioheat Equation is a linear PDE, the reduced-order model is also assumed to be a linear dynamical system. Therefore, a generic reduced-order model for Pennes’ Equation can be written in the form:

    Eq. 11

where the matrix is an inferred linear operator for the ODE system and the matrix is an inferred operator on the input vector to the ODE system. In the context of the examples studied in this article, the vector u will contain a constant term, any time-varying boundary fluxes, and any the time-varying heat sources.

As shown by Peherstorfer and Wilcox [7], the matrices C and D can be obtained by fitting of the dynamics of a full space thermal simulation to the model given by (11). Given the values of m temporal coefficients through time (e.g., those corresponding to the first m spatial modes that were obtained by performing POD on the original simulation data), their time derivatives (estimated by e.g., numerical differencing), and the vector u(t) that was used to perform the simulation, the matrices C and D in (11) can be directly estimated by solving an ordinary least squares problem.

Thermal Dose Modeling with ROMs

The goal of the modeling in this article is to solve for the thermal dose received from a powered wearable device by the skin, which requires solution of the PDEs (3) and (4) in the device and the skin. The projection of the governing PDEs onto a reduced space spanned by a set of POD modes yields a simple ODE system that can be solved more efficiently than the full model. The suggested workflow to solve (3) and (4) and calculate (and validate) the thermal dose received by the skin is as follows:

1. Perform a representative training simulation using a full-space model to obtain T (t, x, y, z) data.
2. Perform POD (e.g., using the Method of Snapshots [6]) on the temperature data to extract the spatial modes.
3. Choose the dimensionality m of the reduced space (e.g., according to the criteria (8) and retain the m spatial modes with the largest eigenvalues.
4. Compute time-independent terms in the matrices and vectors of either (9) or (11) in the m-dimensional reduced space.
5. Choose new dynamics for time-varying source and boundary terms and a new simulation duration.
6. Integrate (9) or (11) using an ODE solver, such as in Matlab, updating time-dependent source and boundary terms as needed at each timestep.
7. Reconstruct the temperature field from the spatial modes and the temporal coefficients using (5).
8. (For validation) Perform an additional full-space simulation to verify the accuracy of the ROM temperature field.

Once the long duration simulation is complete, the CEM43°C value for the exposure can be calculated from (1) and compared to threshold values for the skin tissues.

Illustrative Example

This example considers a long-duration (8-hour) one-dimensional transient simulation of a powered device in contact with the skin. The device and skin are modeled as a multilayered medium with the materials and their respective thicknesses given in Table 1, listed in order from outermost to innermost layer. In the unmodeled spatial dimensions, the device is assumed to be 2.5 cm by 2.5 cm in size.

Table 1: Device and skin layers and their thicknesses

There are two heat generating components within the device: the chip and the battery. For the purposes of this illustrative example, the heat generated (in Watts) by these components over eight hours of continuous use were given by (12) and (13):

    Eq. 12

    Eq. 13

These power profiles are plotted in Figure 2, expressed as an equivalent heat flux based on the assumed cross-sectional area (6.25 × 10⁻⁴m²) of the device.

The outermost glass layer of the device was assumed to be in contact with ambient temperature air at all times. The boundary condition at this surface was a convective heat transfer with a constant heat transfer coefficient of 5 W/m²-K and a time-varying ambient temperature profile, T_amb, given by:

    Eq. 14

The boundary deep within the skin at the end of the modeled inner tissue was represented by constant temperature conditions at the temperature of the blood within the perfuse tissue. The temperature of blood is assumed to be constant at 37°C. The physical properties for the materials and tissues used in this Example are given in Table 2.

Table 2: Thermophysical properties for the materials and tissues in the example

All simulations were performed using code written in Matlab R2018b on a Windows laptop with an i7-7600 CPU (2.8 GHz) and 16 GB of RAM.

Following the analysis procedure outlined in the previous section, the first step was to perform a training simulation to generate representative temperature data. In this case, a 15-minute simulation was performed using greatly simplified boundary conditions and heat generation profiles. Therefore, instead of the power profiles given in (12) and (13), constant heat generation rates of 0.2 W and 0.6 W were assigned to the chip and battery, respectively. Similarly, a constant ambient temperature of 30°C was modeled instead of the profile given in (14). This training simulation was performed using a finite-volume implementation of the governing equations on a grid of 2,750 elements of equal size (10⁻⁵ m) and a constant timestep of 0.5 seconds.

Next, the method of snapshots was used to extract the POD modes. The term U₀ in (5) was chosen to be a constant vector with value 37°C at all points in space to homogenize the inner tissue boundary. The fraction of the total variance captured by each of the first 50 modes is shown in Figure 3. As expected, there is rapid decay after the first few POD modes.

Choosing ε_tol = 10⁻⁶ in (8), the reduced space is formed using the first four POD modes. These spatial modes are plotted in Figure 4.

Intrusive Projection

An intrusive method of projection was considered first and the matrices in (9) were computed using the Galerkin method. The resulting system of four ODEs was then solved using in Matlab for the full eight-hour duration using the power profiles given in (12) and (13) and the ambient temperature profile given in (14). The solution of the ODE system over this entire time horizon takes approximately 1.5 seconds on the hardware indicated previously. The resulting temperature profiles for the surface of the device that is in contact with the epidermis and the basal layer are shown in Figure 5.

The full simulation was then performed using the finite volume method for validation. Resolving the full eight-hour profile with the full order method takes approximately 15 seconds, an order of magnitude longer. The CEM43°C profile was calculated using (1) from both the reduced order model temperature field and the full order model temperature field. The profiles are compared in Figure 6. As Figure 6 shows, the two simulations are in excellent agreement. After eight hours of exposure, the predicted CEM43°C values between the two methods differ by only 0.014 minutes (0.04% relative error for the reduced order method).

This simulation result also shows that even though the surface temperature of the device in contact with the skin exceeded 43°C for over 15 minutes, the total thermal dose as measured by CEM43°C is only 38.4 minutes after eight hours, which is far below the threshold for injurious conditions.

Non-intrusive Projection

The thermal dose analysis was also performed using the non-intrusive projection method. When using this method, some additional requirements must be placed on the training simulation to ensure that the dynamics of the model can be learned from the data. In the intrusive method considered previously, the training simulation was performed with a constant ambient temperature and constant power generation in the heat sources. With the non-intrusive method, simple dynamics were included in the training simulation to differentiate the impacts of the two heat sources, the convective boundary condition, and the time-invariant terms in the model.

For this example, the total duration of the training simulation in the non-intrusive case was also increased to 30 minutes. The chip power was set to increase linearly from 0.1 W to 0.2 W over 15 minutes, and then decrease back to 0.1 W over 15 minutes. Similarly, the battery power increased linearly from 0.3 W to 0.6 W over 15 minutes, and then decreased back to 0.3 W over 15 minutes. The ambient temperature varied quadratically between 27.5°C and 32.5°C according to T_train(t) = 27.5 + 20 ·  ((t − 900)/1800)².
Note that these dynamics are still substantially simpler than those in the full duration simulation.

As with the previous training simulation, the POD modes were calculated from the temperature data and a reduced order model was built using the first four spatial modes. The ODE system was constructed by fitting the training simulation data to (11) using linear least squares. The reduced-order model for the eight-hour simulation was then solved as before using , which took approximately 1.1 seconds on the indicated hardware. A comparison of the predicted CEM43°C values from this method and the finite volume method is shown in Figure 7. At the end of the eight-hour simulation, the difference between the predicted CEM43°C value between the two methods was 0.51 minutes (1.3% relative error for the reduced order method).

While the non-intrusive method required a slightly more complex training simulation and gave a marginally less accurate result, it must be emphasized that the construction of the ODE system using this method is trivial. There is no need to perform any numerical integration to construct the matrices for the dynamical system (11). No information about the geometry and physical properties of modeled system is needed after the training simulation data has been obtained.

Future Work

While the examples shown in this article have focused on a one-dimensional simulation, both the intrusive and non-intrusive methods generalize directly to two and three spatial dimensions. The potential reduction in simulation time is substantial in higher dimensions because the size of the reduced-order model is controlled by the number of modes retained and not the mesh size of the domain. However, practical considerations such as simulation domain size, data storage, and effective training simulations become more challenging in higher dimensions. In future studies, higher dimensional simulations will be performed following the analysis technique established in this article. The implementation of reduced-order models in device control algorithms to understand and forecast thermal shutdown or power throttling requirements for user safety is also an area for future study and development.

Conclusions

In this article, methods have been developed for evaluating thermal doses received by a user from a powered wearable device. These methods do not require repeatedly finding the solution of a complex PDE model to account for changing boundary conditions or heat sources. Once a ROM has been generated using these methods, it can be solved much more efficiently that the full-space PDE model. The predictions from these ROMs can be highly accurate even when evaluating thermal responses to significantly different or more complex dynamics than those on which they were trained. The solution of this heat transfer problem can therefore be divided into two steps: an “offline” step in which the ROM is first generated, and then an “online” step in which the ROM can be run repeatedly, accurately, and efficiently.

References

1. International Data Corporation, Wearables Deliver Double-Digit Growth for Both Q4 and the Full Year 2021, According to IDC, March 8, 2022.
2. ASTM International, C1055-03: Standard Guide for Heated System Surface Conditions that Produce Contact Burn Injuries, 2003.
3. S. A. Sapareto and W.C. Dewey, “Thermal dose determination in cancer therapy,” International Journal of Radiation Oncology, Biology, Physics. 10:787–800, 1984.
4. A.R. Mortiz and F.C. Henriques, “Studies of Thermal Injury: II. The Relative Importance of Time and Surface Temperature in the Causation of Cutaneous Burns,” The American Journal of Pathology, 23(5):695-720, 1947.
5. H. H. Pennes, “Analysis of Tissue and Arterial Temperature in the Resting Human Forearm,” Journal of Applied Physiology, 1:93-122, 1948.
6. L. Sirovich, “Turbulence and the dynamics of coherent structures. I. coherent structures” Quarterly of Applied Mathematics, 45(3):561-571, 1987.
7. B. Peherstorfer and K. Willcox, “Data-driven operator inference for nonintrusive projection-based model reduction,” Computer Methods in Applied Mechanics and Engineering, 306:196-215, 2016.