Hongyun Wang^1*, Parthiv Seetharaman², Shannon E. Foley³, Hong Zhou^4
1Department of Applied Mathematics, University of California, Santa Cruz, CA, USA.
²Department of Mathematics, California Polytechnic State University, San Luis Obispo, CA, USA.
³U.S. Department of Defense, Joint Intermediate Force Capabilities Office, Quantico, VA, USA.
⁴Department of Applied Mathematics, Naval Postgraduate School, Monterey, CA, USA.
DOI: 10.4236/jamp.2026.148155   PDF    HTML   XML   8 Downloads   54 Views   Citations

Abstract

In electromagnetic heating and other applications, it is necessary to compute the volume enclosed by an isosurface of the 3D temperature distribution that is numerically represented on a rectangular grid. This situation arises naturally when the temperature distribution is obtained by solving a partial differential equation numerically using a finite difference method. Given the temperature distribution on the 3D grid, the isosurface of a prescribed value is represented approximately by a triangulation, a collection of triangles with vertices on the grid lines. The vertices are determined by a linear interpolation to approximate the locations where the temperature is at the prescribed level. The region enclosed by the isosurface is approximated by that enclosed by the triangulation, which is a set of tetrahedra. The enclosed volume is approximated by summing the contribution of these tetrahedra. This volume approximation is analogous to approximating a curve using line segments and is limited to second-order accuracy. In this study, we combine extrapolation with the triangulation approximation to develop a more accurate method for computing the volume of the region in which the temperature is above a given threshold.

Keywords

Electromagnetic Heating, Skin Activated Volume, Isosurface, Triangulation, Extrapolation

Share and Cite:

1. Introduction

In electromagnetic heating of the human body, the skin surface is exposed to millimeter waves (MMW) in the frequency range of 30 - 300 GHz. The millimeter-wave band corresponds to wavelengths between 1 mm at 300 GHz and 10 mm at 30 GHz. MMWs have a wide range of applications [1]-[7], including 5G networks [8], detection of explosives [9], and autonomous railway systems [10]. Given their increasing prevalence, it is important to study the biological effects of MMW exposure on humans. The primary effect of MMW exposure on humans is electromagnetic heating. Skin is a lossy medium. MMW wave is rapidly attenuated during propagation in skin, and the electromagnetic power is absorbed within a thin layer of skin, which becomes a heat source in the skin thermal evolution. When the local temperature reaches the activation threshold, thermal nociceptors at that skin location are activated, inducing a heating sensation. Over an extended duration of MMW exposure, the skin temperature keeps increasing, which activates nociceptors over a larger region and intensifies the induced heating sensation. For the purpose of assessing the induced heating sensation and/or assessing the thermal injury risk, it is necessary to compute the skin volume where the local temperature is at or above a given value. This skin volume corresponds to the region enclosed by the isosurface of the 3D skin temperature at the given level. In this study, we cast this task into a broad mathematical setting and consider the general problem of computing the 3D volume enclosed by an isosurface of function $f\left(x,y,z\right)$ . We focus on the practical situation where function $f\left(x,y,z\right)$ is described only discretely on a 3D rectangular grid. This situation arises naturally in many thermal applications where the temperature distribution is obtained by solving the heat equation numerically on a 3D grid.

The remainder of the paper is organized as follows. Section 2 presents the mathematical formulation of the method for computing the enclosed volume and for carrying out extrapolation. Section 3 examines the accuracy of the method, respectively, without and with extrapolation, on several test problems in which the exact enclosed volume is known analytically. Section 4 tests the performance of the methods on a simplified case of electromagnetic heating where an idealized skin tissue of uniform material properties is exposed to an axial-symmetric Gaussian beam at perpendicular incidence, and the lateral heat conduction is neglected. In this idealized electromagnetic heating problem, the activated volume can be practically calculated up to machine precision by utilizing the particular axial symmetry in this model problem. Section 5 summarizes the main results of the study.

2. Mathematical Formulation

Consider region $D$ in ${ℝ}^3$ . We view it as enclosed by its boundary surface $S=\partial D$ . Let $x_0=\left(x_0,y_0,z_0\right)$ be a point in ${ℝ}^3$ . We use the divergence theorem and the fact that $\nabla \cdot \left(x-x_0\right)=3$ to connect a surface integral to the volume of region $D$ .

$\underset{\text{divergence theorem}}{\underbrace{{\displaystyle {\oint }_{\partial D}\left(x-x_0\right)\cdot \text{d}S}={\displaystyle {\int }_D\nabla \cdot \left(x-x_0\right)\text{d}V}}}=3\times \text{volume}\left(D\right)$ (1)

where $\text{d}S=\widehat{n}\text{d}S$ is the vector area element, containing both the magnitude of the area element $\text{d}S$ and its outward unit normal direction $\widehat{n}$ . From (1), we express the volume of region $D$ in terms of a surface integral over its boundary surface.

$\text{volume}\left(D\right)=\frac{1}{3}{\displaystyle {\oint }_S\left(x-x_0\right)\cdot \text{d}S},x_0\in {ℝ}^3$ (2)

where $S=\partial D$ is the enclosing boundary surface. (2) is valid no matter whether $x_0$ is inside or outside region $D$ or on its boundary surface.

In electromagnetic heating applications, region $D$ is defined as having temperature at or above a given level. As a result, its boundary surface $S=\partial D$ is the isosurface of the 3D skin temperature at the given level. Let $f\left(x\right)$ be the skin temperature at position $x$ . We consider the situation where function $f\left(x\right)$ is only represented by its values on a 3D rectangular grid. In standard computational packages (such as MATLAB), given a discrete representation of function $f\left(x\right)$ and a level value, the isosurface $S$ is approximated by a triangulation representation, a collection of triangular patches, as illustrated in Figure 1.

$S={\displaystyle \sum }_{j=1}^N{S}_j\approx {\displaystyle \sum }_{j=1}^N\underset{\text{planar patch approximation of}{S}_j}{\underbrace{\text{Δ}\left(x_1^{\left(j\right)},x_2^{\left(j\right)},x_3^{\left(j\right)}\right)}}$ (3)

In (3), $S_j$ represents a small piece of the true curly surface $S$ while $\text{Δ}\left(x_1^{\left(j\right)},x_2^{\left(j\right)},x_3^{\left(j\right)}\right)$ is a triangular patch approximating $S_j$ . $N$ is the number of patches in the triangulation. The surface integral over $S$ is approximated by the sum of integrals over individual patches.

${\displaystyle {\int }_S\left(x-x_0\right)\cdot \text{d}S}\approx {\displaystyle \sum _{j=1}^N{\displaystyle {\int }_{\Delta \left(x_1^{\left(j\right)},x_2^{\left(j\right)},x_3^{\left(j\right)}\right)}\left(x-x_0\right)\cdot \text{d}S}}$ (4)

Each patch is a planar triangle. The integral over a patch has an analytical expression

${\displaystyle {\int }_{\Delta \left(x_1^{\left(j\right)},x_2^{\left(j\right)},x_3^{\left(j\right)}\right)}\left(x-x_0\right)\cdot \text{d}S}=\frac{1}{2}\left(\left(x_1^{\left(j\right)}-x_0\right)\times \left(x_2^{\left(j\right)}-x_0\right)\right)\cdot \left(x_3^{\left(j\right)}-x_0\right)$ (5)

Therefore, the volume of region $D$ (the region enclosed by surface $S$ ) is approximated by

(a) (b)

Figure 1. Top: region enclosed by the triangulation of an isosurface. Bottom: detailed view of a triangle element in the isosurface triangulation.

$\boxed{\text{volume}\left(D\right)\approx \underset{\text{volume enclosed by the triangulation}}{\underbrace{\frac{1}{6}{\displaystyle \sum }_{j=1}^N\left(\left(x_1^{\left(j\right)}-x_0\right)\times \left(x_2^{\left(j\right)}-x_0\right)\right)\cdot \left(x_3^{\left(j\right)}-x_0\right)}}\equiv V\left(\Delta x\right)}$ (6)

where $V\left(\text{Δ}x\right)$ denotes the volume enclosed by the triangulation of the isosurface based on the discrete representation of function $f\left(x\right)$ on the grid of mesh spacing $\text{Δ}x$ . $V\left(\text{Δ}x\right)$ is used as an approximation to the volume enclosed by the true isosurface, which is not known exactly and has some uncertainty if function $f\left(x\right)$ is represented only on a discrete grid. In particular, the notation $V\left(\text{Δ}x\right)$ indicates the effect of mesh spacing $\text{Δ}x$ on the isosurface triangulation and on the volume approximation explicitly.

In the triangulation, the true isosurface in the 3D is replaced with a set of planar triangle elements. This is analogous to the situation where a curve in 3D is replaced with a set of line segments, which yields a second-order approximation. Thus, we expect $V\left(\text{Δ}x\right)$ to be a second-order approximation of the true volume $V_{\text{true}}\equiv \text{volume}\left(D\right)$ .

$V\left(\Delta x\right)=V_{\text{true}}+C_2{\left(\Delta x\right)}^2+o\left({\left(\Delta x\right)}^2\right)$ (7)

Given the discrete representation of function $f\left(x\right)$ on a 3D grid of mesh spacing $\text{Δ}x$ , we can always down-sample to obtain a coarse representation of function $f\left(x\right)$ on a grid of mesh spacing $\left(2\text{Δ}x\right)$ . From the coarse representation, we construct a coarse triangulation of the isosurface and the corresponding volume approximation $V\left(2\text{Δ}x\right)$ , which satisfies

$V\left(2\Delta x\right)=V_{\text{true}}+C_2{\left(2\Delta x\right)}^2+o\left({\left(\Delta x\right)}^2\right)$ (8)

In (7), the dominant part of error in $V\left(\text{Δ}x\right)$ is $C_2{\left(\text{Δ}x\right)}^2$ . (8) based on the coarse representation, it contains a similar term $C_2{\left(2\text{Δ}x\right)}^2$ . We combine (7) and (8) to eliminate the dominant error $C_2{\left(\text{Δ}x\right)}^2$ in $V\left(\text{Δ}x\right)$ . The result is a more accurate approximation to $V_{\text{true}}$ .

$\boxed{V_{\text{extrap}}\left(\Delta x\right)\equiv \frac{4V\left(\Delta x\right)-V\left(2\Delta x\right)}{3}=V_{\text{true}}+o\left({\left(\Delta x\right)}^2\right)}$ (9)

Mathematically, this is an extrapolation: we are using $V\left(2\text{Δ}x\right)$ and $V\left(\text{Δ}x\right)$ to predict $V\left(0_{+}\right)$ , the numerical approximation at an infinitesimal mesh spacing.

3. Several Test Problems of Enclosed Volumes with Analytical Solutions

We study the numerical accuracy of the volume computation method (6) and the extrapolation improvement (9). To examine the true numerical error, we select several test problems in which the volume enclosed by the given isosurface has a closed-form analytical expression. In the discussion below, let $f\left(x\right)$ be the underlying function, mimicking the role of the 3D skin temperature distribution in the electromagnetic heating problem. Let $D$ be the region enclosed by the isosurface of function $f\left(x\right)$ at a given value $f_c$ . As introduced in the last section, $\text{volume}\left(D\right)$ denotes the true volume of the enclosed region $D$ , and $V\left(\text{Δ}x\right)$ and $V_{\text{extrap}}\left(\text{Δ}x\right)$ denote the numerical volumes of respectively method (6) and method (9) based on a discrete representation of $f\left(x\right)$ on a grid of mesh spacing $\text{Δ}x$ . Method (9) is built on method (6) with extrapolation incorporated.

The true numerical errors, respectively, in $V\left(\text{Δ}x\right)$ and in $V_{\text{extrap}}\left(\text{Δ}x\right)$ are defined as

$\{\begin{array}{l}\left(\text{error}\text{in}V\left(\Delta x\right)\right)=\left|V\left(\Delta x\right)-\text{volume}\left(D\right)\right|\\ \left(\text{error}\text{in}{V}_{\text{extrap}}\left(\Delta x\right)\right)=\left|V_{\text{extrap}}\left(\Delta x\right)-\text{volume}\left(D\right)\right|\end{array}$ (10)

3.1. The Unit Sphere

$\{\begin{array}{l}\text{function}:f\left(x\right)=1-\left(x^2+y^2+z^2\right),x=\left(x,y,z\right)\\ \text{enclosed region}:D=\left\{\left(x,y,z\right)|f\left(x\right)\ge 0\right\}\\ \text{volume}\left(D\right)=\frac{4\pi }{3}\end{array}$ (11)

Numerical results for the unit sphere are shown in Figure 2. As predicted in (7), the error in $V\left(\text{Δ}x\right)$ is approximately proportional to ${\left(\text{Δ}x\right)}^2$ . The error in $V_{\text{extrap}}\left(\text{Δ}x\right)$ is significantly below that in $V\left(\text{Δ}x\right)$ , demonstrating the benefit of extrapolation. Note that the extrapolation component does not require additional data or information on function $f\left(x\right)$ . $V\left(2\text{Δ}x\right)$ in the extrapolation is obtained from down-sampling the given discrete representation of $f\left(x\right)$ .

(a) (b)

Figure 2. Left: Triangulation of the unit sphere as an isosurface of function $f\left(x\right)$ given in (11). Right: errors in numerical volumes $V\left(\text{Δ}x\right)$ and $V_{\text{extrap}}\left(\text{Δ}x\right)$ vs mesh spacing $\text{Δ}x$ . Dashed line is ${\left(\text{Δ}x\right)}^2$ vs $\text{Δ}x$ serving as a visual guide.

3.2. An Axis-Aligned Ellipsoid

$\{\begin{array}{l}\text{Function}:f\left(x\right)=1-\left(\frac{x^2}{a^2}+\frac{y^2}{b^2}+\frac{z^2}{c^2}\right),x=\left(x,y,z\right)\\ \text{where}a=1,b=\frac{1}{4},c=\frac{1}{2}\\ \text{enclosed region}:D=\left\{\left(x,y,z\right)|f\left(x\right)\ge 0\right\}\\ \text{volume}\left(D\right)=\frac{4\pi }{3}abc\end{array}$ (12)

In Figure 3 for the axis-aligned ellipsoid, the volume computation methods show the same behavior as in Figure 2. The error in $V\left(\text{Δ}x\right)$ is approximately proportional to ${\left(\text{Δ}x\right)}^2$ ; the error in $V_{\text{extrap}}\left(\text{Δ}x\right)$ is significantly smaller than that in $V\left(\text{Δ}x\right)$ , confirming the second-order accuracy of triangulation and the accuracy improvement of extrapolation.

(a) (b)

Figure 3. Left: Triangulation of the axis-aligned ellipsoid as an isosurface of $f\left(x\right)$ in (12). Right: errors in numerical volumes $V\left(\text{Δ}x\right)$ and $V_{\text{extrap}}\left(\text{Δ}x\right)$ vs mesh spacing $\text{Δ}x$ .

3.3. A Rotated Ellipsoid

$\{\begin{array}{l}\text{Function}:f\left(x\right)=1-\left(\frac{{\tilde{x}}^2}{a^2}+\frac{{\tilde{y}}^2}{b^2}+\frac{{\tilde{z}}^2}{c^2}\right),x=\left(x,y,z\right)\\ \text{where}\left(\begin{array}{l}\tilde{x}\hfill \\ \tilde{y}\hfill \\ \tilde{z}\hfill \end{array}\right)=A\left(\begin{array}{l}x\hfill \\ y\hfill \\ z\hfill \end{array}\right),A=\left(\begin{array}{ccc}\frac{\sqrt{3}}{2\sqrt{2}}& \frac{1}{\sqrt{2}}& \frac{1}{2\sqrt{2}}\\ \frac{-\sqrt{3}}{2\sqrt{2}}& \frac{1}{\sqrt{2}}& \frac{-1}{2\sqrt{2}}\\ \frac{-1}{2}& 0& \frac{\sqrt{3}}{2}\end{array}\right)\\ a=1,b=\frac{1}{4},c=\frac{1}{2}\\ \text{enclosed region}:D=\left\{\left(x,y,z\right)|f\left(x\right)\ge 0\right\}\\ \text{volume}\left(D\right)=\frac{4\pi }{3}abc\end{array}$(13)

In (13), matrix $A$ is orthogonal, satisfying $A^{\text{T}}A=I$ . Geometrically, multiplying by $A$ corresponds to rotating the 3D region with respect to the z-axis by 45˚ and then rotating it with respect to the y-axis by (−30˚). The enclosed region in (13) is congruent to that of Subsection 3.2 via an orthogonal transformation, which preserves the analytical volume.

The discrete representation of function $f\left(x\right)$ on a fixed 3D grid, however, varies with the orientation. Accordingly, the corresponding numerical volume based on the discrete representation is expected to be affected. In this sense, the rotated ellipsoid provides a completely new test on the volume computation methods, although the exact volume stays the same. For the rotated ellipsoid, in Figure 4, we observe the same behavior as in Figure 2 and Figure 3. The error in $V\left(\text{Δ}x\right)$ is approximately proportional to ${\left(\text{Δ}x\right)}^2$ ; the error in $V_{\text{extrap}}\left(\text{Δ}x\right)$ is much lower than that in $V\left(\text{Δ}x\right)$ .

(a) (b)

Figure 4. Left: Triangulation of the rotated ellipsoid as an isosurface of $f\left(x\right)$ in (13). Right: errors in numerical volumes $V\left(\text{Δ}x\right)$ and $V_{\text{extrap}}\left(\text{Δ}x\right)$ vs mesh spacing $\text{Δ}x$ .

3.4. An Axis-Aligned Superellipsoid

$\{\begin{array}{l}\text{Function}:f\left(x\right)=1-\left({\left[{\left(\frac{x}{a}\right)}^{\frac{2}{{\epsilon }_2}}+{\left(\frac{y}{b}\right)}^{\frac{2}{{\epsilon }_2}}\right]}^{\frac{{\epsilon }_2}{{\epsilon }_1}}+{\left(\frac{z}{c}\right)}^{\frac{2}{{\epsilon }_1}}\right),x=\left(x,y,z\right)\\ \text{where}{\epsilon }_1={\epsilon }_2=\frac{3}{4},a=1,b=\frac{1}{4},c=\frac{1}{2}\\ \text{enclosed region}:D=\left\{\left(x,y,z\right)|f\left(x\right)\ge 0\right\}\\ \text{volume}\left(D\right)=2abc{\epsilon }_1{\epsilon }_2\beta \left(\frac{{\epsilon }_1}{2},{\epsilon }_1+1\right)\beta \left(\frac{{\epsilon }_2}{2},{\epsilon }_2+2\right)\end{array}$ (14)

where $\beta \left(a,b\right)$ is the beta function defined as

$\beta \left(a,b\right)\equiv {\displaystyle {\int }_0^1{t}^{a-1}{\left(1-t\right)}^{b-1}\text{d}t}$

The beta function is available in standard software packages.

The regular ellipsoid (12) is a special case of the superellipsoid (14) at ${\epsilon }_1={\epsilon }_2=1$ . In the case of ${\epsilon }_1={\epsilon }_2=\frac{2}{n}$ , the enclosed region is described by

${\left(\frac{x}{a}\right)}^n+{\left(\frac{y}{b}\right)}^n+{\left(\frac{z}{c}\right)}^n\le 1$

As $n\to +\infty$ , the enclosed region converges to the 3D rectangular prism:

$\left[-a,a\right]\times \left[-b,b\right]\times \left[-c,c\right]$

For ${\epsilon }_1={\epsilon }_2<1$ , the enclosed superellipsoid is between the regular ellipsoid and the 3D rectangular prism. Indeed, the principal sides of the superellipsoid shown in the left panel of Figure 5 are flatter than those of the regular ellipsoid in Figure 3. For the superellipsoid, even with this shape change in the enclosed region, the error trends in volume computation remain the same as in Figures 2-4. The error in $V\left(\text{Δ}x\right)$ is proportional to ${\left(\text{Δ}x\right)}^2$ ; the error in $V_{\text{extrap}}\left(\text{Δ}x\right)$ is significantly smaller than that in $V\left(\text{Δ}x\right)$ .

(a) (b)

Figure 5. Left: Triangulation of the axis-aligned superellipsoid as an isosurface of $f\left(x\right)$ in (14). Right: errors in numerical volumes $V\left(\text{Δ}x\right)$ and $V_{\text{extrap}}\left(\text{Δ}x\right)$ vs mesh spacing $\text{Δ}x$ .

3.5. A Rotated Superellipsoid

$\{\begin{array}{l}\text{Function}:f\left(x\right)=1-\left({\left[{\left(\frac{\tilde{x}}{a}\right)}^{\frac{2}{{\epsilon }_2}}+{\left(\frac{\tilde{y}}{b}\right)}^{\frac{2}{{\epsilon }_2}}\right]}^{\frac{{\epsilon }_2}{{\epsilon }_1}}+{\left(\frac{\tilde{z}}{c}\right)}^{\frac{2}{{\epsilon }_1}}\right),x=\left(x,y,z\right)\\ \text{where}\left(\begin{array}{l}\tilde{x}\hfill \\ \tilde{y}\hfill \\ \tilde{z}\hfill \end{array}\right)=A\left(\begin{array}{l}x\hfill \\ y\hfill \\ z\hfill \end{array}\right),A=\left(\begin{array}{ccc}\frac{\sqrt{3}}{2\sqrt{2}}& \frac{1}{\sqrt{2}}& \frac{1}{2\sqrt{2}}\\ \frac{-\sqrt{3}}{2\sqrt{2}}& \frac{1}{\sqrt{2}}& \frac{-1}{2\sqrt{2}}\\ \frac{-1}{2}& 0& \frac{\sqrt{3}}{2}\end{array}\right)\\ {\epsilon }_1={\epsilon }_2=\frac{3}{4},a=1,b=\frac{1}{4},c=\frac{1}{2}\\ \text{enclosed region}:D=\left\{\left(x,y,z\right)|f\left(x\right)\ge 0\right\}\\ \text{volume}\left(D\right)=2abc{\epsilon }_1{\epsilon }_2\beta \left(\frac{{\epsilon }_1}{2},{\epsilon }_1+1\right)\beta \left(\frac{{\epsilon }_2}{2},{\epsilon }_2+2\right)\end{array}$(15)

The enclosed region in (15) is congruent to that in (14) of Subsection 3.4 via an orthogonal transformation. In (15), the orthogonal transformation represented by matrix $A$ is the same as that in (13). Similar to what we observed in Figures 2-5, Figure 6 shows that the errors in $V\left(\text{Δ}x\right)$ is proportional to ${\left(\text{Δ}x\right)}^2$ ; the error in $V_{\text{extrap}}\left(\text{Δ}x\right)$ is significantly smaller than that in $V\left(\text{Δ}x\right)$ , demonstrating the steady improvement of extrapolation over several test problems.

(a) (b)

Figure 6. Left: Triangulation of the rotated superellipsoid as an isosurface of $f\left(x\right)$ in (15). Right: errors in numerical volumes $V\left(\text{Δ}x\right)$ and $V_{\text{extrap}}\left(\text{Δ}x\right)$ vs mesh spacing $\text{Δ}x$ .

3.6. An Axis-Aligned Elliptic Ring Torus

$\{\begin{array}{l}\text{function}:f\left(x\right)=r^2-\left({\left[\sqrt{x^2+y^2}-R\right]}^2+{\left(\frac{z}{c}\right)}^2\right),x=\left(x,y,z\right)\\ \text{where}R=0.8,r=0.28,c=\frac{1}{0.7}\\ \text{enclosed region}:D=\left\{\left(x,y,z\right)|f\left(x\right)\ge 0\right\}\\ \text{volume}\left(D\right)=2{\pi }^{2}c{r}^{2}R\end{array}$ (16)

A triangulation of the axis-aligned elliptic ring torus is shown in the left panel of Figure 7. Geometrically, the centerline of the torus is the circle of radius $R$ centered at the origin in the $\left(x,y\right)$ -plane. Perpendicular to the centerline, the cross section of the torus is an ellipse with the semi minor axis $r$ in the horizontal direction and the semi major axis ($cr$ ) in the z-direction. For the elliptic ring torus in (16), $c>1$ . The exact volume of the enclosed region is given by the product of the center line arclength ($2\pi R$ ) and the cross-sectional area ($\pi c{r}^2$ ). We use the ring torus as a test problem to examine the error trends in numerical volumes.

(a) (b)

Figure 7. Left: Triangulation of the axis-aligned elliptic ring torus as an isosurface of $f\left(x\right)$ in (16). Right: errors in numerical volumes $V\left(\text{Δ}x\right)$ and $V_{\text{extrap}}\left(\text{Δ}x\right)$ vs mesh spacing $\text{Δ}x$ .

Note that numerical volumes are computed from a discrete representation of function $f\left(x\right)$ on a 3D grid of mesh spacing $\text{Δ}x$ , not using the full analytical expression of $f\left(x\right)$ . The exact volume is not used in computing numerical volumes; it is used only in evaluating the errors of numerical volumes. The right panel of Figure 7 shows that the errors in $V\left(\text{Δ}x\right)$ and $V_{\text{extrap}}\left(\text{Δ}x\right)$ both decrease with $\text{Δ}x$ with the latter significantly below the former.

3.7. A Rotated Elliptic Ring Torus

$\{\begin{array}{l}\text{function}:f\left(x\right)=r^2-\left({\left[\sqrt{{\tilde{x}}^2+{\tilde{y}}^2}-R\right]}^2+{\left(\frac{\tilde{z}}{c}\right)}^2\right),x=\left(x,y,z\right)\\ \text{where}\left(\begin{array}{l}\tilde{x}\hfill \\ \tilde{y}\hfill \\ \tilde{z}\hfill \end{array}\right)=B\left(\begin{array}{l}x\hfill \\ y\hfill \\ z\hfill \end{array}\right),B=\left(\begin{array}{ccc}\frac{\sqrt{3}}{2\sqrt{2}}& \frac{\sqrt{3}}{2\sqrt{2}}& \frac{1}{2}\\ \frac{-1}{\sqrt{2}}& \frac{1}{\sqrt{2}}& 0\\ \frac{-1}{2\sqrt{2}}& \frac{-1}{2\sqrt{2}}& \frac{\sqrt{3}}{2}\end{array}\right)\\ R=0.8,r=0.28,c=\frac{1}{0.7}\\ \text{enclosed region}:D=\left\{\left(x,y,z\right)|f\left(x\right)\ge 0\right\}\\ \text{volume}\left(D\right)=2{\pi }^{2}c{r}^{2}R\end{array}$(17)

In (17), the orthogonal transformation represented by matrix $B$ is different from that by matrix A in (13) and (15). If we first rotate the elliptic ring torus about the z-axis, due to the axial symmetry of the object, this rotation would not change the object at all. Instead, in (17), the transformation represented by matrix $B$ corresponds to rotating the 3D object first with respect to the y-axis by (−30˚) and then rotating it with respect to the z-axis by 45˚. A triangulation of the rotated elliptic ring torus is shown in the left panel of Figure 8. The error trends of numerical volumes are plotted in the right panel of Figure 8. Similar to what we observed in all other test problems, the error in $V\left(\text{Δ}x\right)$ is proportional to ${\left(\text{Δ}x\right)}^2$ ; the error in $V_{\text{extrap}}\left(\text{Δ}x\right)$ is significantly below that in $V\left(\text{Δ}x\right)$ .

4. A Simplified Case of Electromagnetic Heating

We study the electromagnetic heating of a flat skin by an incident electromagnetic wave of millimeter wavelength. We consider the idealized situation where the incident wave is axially symmetric and is perpendicular to the flat skin surface, as illustrated in Figure 9. We first establish the coordinate system as described in our previous study [11] [12]. The normal direction of the flat skin surface pointing into the skin tissue is selected as the positive z-direction. The skin surface is at $z=0$ , and $z>0$ corresponds to the skin tissue beneath the surface. The center of the incident beam is selected as the origin of the $\left(x,y\right)$ -plane. In the skin tissue, the 3D coordinates of a point are represented by $\left(x,y,z\right)$ , $z\ge 0$ .

(a) (b)

Figure 8. Left: Triangulation of the rotated elliptic ring torus as an isosurface of $f\left(x\right)$ in (17). Right: errors in numerical volumes $V\left(\text{Δ}x\right)$ and $V_{\text{extrap}}\left(\text{Δ}x\right)$ vs mesh spacing $\text{Δ}x$ .

Figure 9. The coordinate system in skin electromagnetic heating in the case of flat skin with an incident beam perpendicular to the skin.

4.1. Physical Quantities, Equations, and Non-Dimensionalization

We introduce the physical quantities and equations used in our model.

• $T\left(x,y,z,t\right)$ is the skin temperature at position $\left(x,y,z\right)$ at time $t$ .

• ${\rho }_m$ is the mass density, $C_p$ the specific heat capacity, and $k$ the heat conductivity of skin. In this study, we consider a single layer of uniform skin where all material properties are independent of position $\left(x,y,z\right)$ .

• The energy balance gives the governing equation for $T\left(x,y,z,t\right)$ [13]:

$\underset{\begin{array}{c}\text{rate of change}\\ \text{of heat in skin}\end{array}}{\underbrace{{\rho }_m{C}_p\frac{\partial T}{\partial t}}}=\underset{\begin{array}{c}\text{net heat influx}\\ \text{from conduction}\end{array}}{\underbrace{k\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2}\right)}}+\underset{\begin{array}{c}\text{absorbed power}\\ \text{as heat source}\end{array}}{\underbrace{P_d\left(x,y\right)\mu {\text{e}}^{-\mu z}}}$ (18)

• $\mu$ is the absorption coefficient of skin for the electromagnetic frequency as the wave propagates in the skin. $\mu$ describes the fraction of power absorbed per unit depth as the beam propagates in the skin, which is a lossy medium. The remaining power density at depth $z$ is governed by the Beer-Lambert law [14].

$\left(\text{remaining power density at depth }z\right)\propto {\text{e}}^{-\mu z}$ (19)

$\mu$ has the physical dimension of 1/[length]. ($1/\mu$ ) gives the characteristic penetration depth of the electromagnetic frequency, the depth at which the power density is attenuated by a factor of $1/e$ . ($1/\mu$ ) serves as the depth scale in non-dimensionalization.

• $P_d\left(x,y\right)$ is the power density passing through the surface into skin tissue. Since skin is a lossy medium ($\mu >0$ ), as described in Beer-Lambert law in (19), the remaining power density at depth $z$ is virtually zero at large ($\mu z$ ). All of $P_d\left(x,y\right)$ is converted to heat during the electromagnetic propagation in the skin. The heat production rate per volume at depth $z$ is $P_d\left(x,y\right)\mu {\text{e}}^{-\mu z}$ , which is the heat source in the skin temperature evolution governed by the 3D heat equation (18).

• In our simple model, the power density $P_d\left(x,y\right)$ is an axially symmetric Gaussian.

$P_d\left(x,y\right)=P_d^{\left(0\right)}\exp \left(\frac{-2\left(x^2+y^2\right)}{w^2}\right)=P_d^{\left(0\right)}\exp \left(\frac{-\left(x^2+y^2\right)}{2{\sigma }^2}\right)$

where $w$ is the standard beam radius defined as the distance from the beam center to the location where the power density drops to $1/e^2\approx 13.5\text{\%}$ of the peak value; $\sigma$ is the RMS (root mean square) beam radius, which is the standard deviation of the Gaussian distribution. These two are related by $\sigma =w/2$ .

• In the z-direction, the characteristic scale is given by the penetration depth $z_s=1/\mu$ . In the $\left(x,y\right)$ -lateral directions, the characteristic scale is given by the RMS beam radius $r_s=\sigma$ . In the time direction, the characteristic scale is given by the quantity $t_s\equiv \frac{{\rho }_m{C}_p}{k{\mu }^2}$ . With these scales, we carry out non-dimensionalization and recycle the same notations for all non-dimensional functions and variables afterward. The resulting non-dimensional governing equation is:

$\underset{\text{governing equation after non-dimensionalization}}{\underbrace{\frac{\partial T}{\partial t}=\frac{{\partial }^{2}T}{\partial z^2}+{\epsilon }^2\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}\right)+P_d^{\left(0\right)}\exp \left(\frac{-\left(x^2+y^2\right)}{2}\right){\text{e}}^{-z}}}$ (20)

where $\epsilon \equiv \frac{z_s}{r_s}$ is the ratio of the penetration depth to the lateral scale [15].

• In electromagnetic heating applications using millimeter wavelength, the electromagnetic penetration depth is sub-millimeter while the RMS radius of the Gaussian beam is several centimeters [16]. In these applications, the depth-to-lateral-length ratio satisfies $\epsilon \ll 1$ . Thus, the effect of the lateral heat conduction ${\epsilon }^2\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}\right)$ is negligible in comparison with other terms in (20). To the leading order, the governing equation becomes:

$\underset{\text{simplified non-dimensional equation}}{\underbrace{\frac{\partial T}{\partial t}=\frac{{\partial }^{2}T}{\partial z^2}+P_d^{\left(0\right)}\exp \left(\frac{-\left(x^2+y^2\right)}{2}\right){\text{e}}^{-z}}}$ (21)

4.2. Initial Boundary Value Problem and Its Analytical Solution

We continue discussing the simplified system governing the skin temperature evolution.

• Before the start of electromagnetic exposure, the skin is in its thermal homeostasis. The homeostatic temperature varies gradually from the core temperature inside skin tissue to a lower temperature at the skin surface. In our simplified model, however, we neglect this spatial variation in homeostatic temperature and assume that skin baseline temperature is a uniform constant $T_{\text{base}}^{\left(\text{phy}\right)}$ at the beginning of electromagnetic exposure. This assumption is justified for millimeter wave exposure because the heating is concentrated in a very thin skin layer in which the variation of homeostatic temperature is small. This assumption leads to a simple initial condition.

$T^{\left(\text{phy}\right)}\left(x,y,z,0\right)=T_{\text{base}}^{\left(\text{phy}\right)}$

• In non-dimensionalization (20), physical temperature $T^{\left(\text{phy}\right)}$ is shifted by $T_{\text{base}}^{\left(\text{phy}\right)}$ and normalized by the temperature scale $\text{Δ}T\equiv T_{\text{act}}^{\left(\text{phy}\right)}-T_{\text{base}}^{\left(\text{phy}\right)}$ where $T_{\text{act}}^{\left(\text{phy}\right)}$ is the physical temperature threshold for activating skin thermal nociceptors. Given the physical baseline temperature $T_{\text{base}}^{\left(\text{phy}\right)}$ , the temperature scale $\text{Δ}T$ is the temperature increase needed for skin thermal nociceptor activation. After normalization, the non-dimensional versions of baseline temperature and temperature threshold are

$T_{\text{base}}=0,T_{\text{act}}=1$

• In homeostasis, the skin loses heat to the environment at its surface via blackbody radiation and convective cooling. However, the rate of heat loss from the skin surface is low. In a high-intensity, short-duration electromagnetic exposure, the rate of electromagnetic energy entering the skin tissue is much greater than the rate of slow heat loss at the skin surface. Thus, we neglect the slow heat loss at the skin surface and assume no heat exchange with the environment there. This assumption leads to a simple boundary condition.

${\frac{\partial T\left(x,y,z,t\right)}{\partial z}|}_{z=0}=0$

• The governing system for the evolution of skin temperature consists of a simplified non-dimensional equation and initial and boundary conditions.

$\{\begin{array}{l}\frac{\partial T}{\partial t}=\frac{{\partial }^{2}T}{\partial z^2}+P_d^{\left(0\right)}\exp \left(\frac{-\left(x^2+y^2\right)}{2}\right){\text{e}}^{-z}\\ {\frac{\partial T\left(x,y,z,t\right)}{\partial z}|}_{z=0}=0\\ T\left(x,y,z,0\right)=0\end{array}$ (22)

• Note that (22) has three key features: 1) it is linear; 2) $\left(x,y\right)$ are not involve differential variables; and 3) the source term has separate dependences on $\left(x,y\right)$ and on $z$ . These features imply that the skin temperature $T\left(x,y,z,t\right)$ has the form

$T\left(x,y,z,t\right)=\underset{\left(x,y\right)\text{-dependence}}{\underbrace{P_d^{\left(0\right)}\exp \left(\frac{-\left(x^2+y^2\right)}{2}\right)}}\underset{z\text{-dependence}}{\underbrace{U\left(z,t\right)}}$ (23)

where $U\left(z,t\right)$ is a function of $\left(z,t\right)$ only and satisfies

$\{\begin{array}{l}\frac{\partial U}{\partial t}=\frac{{\partial }^{2}U}{\partial z^2}+{\text{e}}^{-z}\\ {\frac{\partial U\left(z,t\right)}{\partial z}|}_{z=0}=0\\ U\left(z,0\right)=0\end{array}$ (24)

• The solution of IBVP (24) has a closed-form analytical expression [17].

$U\left(z,t\right)=\{\begin{array}{l}\frac{{\text{e}}^{-z+t}}{2}\text{erfc}\left(\frac{-z+2t}{\sqrt{4t}}\right)+\frac{{\text{e}}^{z+t}}{2}\text{erfc}\left(\frac{z+2t}{\sqrt{4t}}\right)\\ -{\text{e}}^{-z}-z\text{erfc}\left(\frac{z}{\sqrt{4t}}\right)+\frac{2\sqrt{t}}{\sqrt{\pi }}{\text{e}}^{\frac{-z^2}{4t}},t>0\\ 0,t\le 0\end{array}$ (25)

where $\text{erfc}\left(s\right)$ is the complementary error function.

$\text{erfc}\left(s\right)=\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_s^{+\infty }\text{exp}\left(-t^2\right)\text{d}t}$

4.3. A Test Problem for Skin-Activated Volume

In this subsection, all quantities are non-dimensional. The activated skin region $D$ at time $t$ is enclosed by the isosurface of skin temperature $T\left(x,y,z,t\right)$ at value $T_{\text{act}}=1$ . Using the expression of $T\left(x,y,z,t\right)$ in (23), we write the activated skin region $D$ as

$\{\begin{array}{l}\text{function}:f\left(x\right)=P_d^{\left(0\right)}\exp \left(\frac{-\left(x^2+y^2\right)}{2}\right)U\left(z,t\right)-1,x=\left(x,y,z\right)\\ \text{enclosed region}:D=\left\{\left(x,y,z\right)|f\left(x\right)\ge 0\right\}\end{array}$ (26)

In this simplified electromagnetic heating problem, the activated region $D$ is axially symmetric. In more realistic electromagnetic heating applications, the electromagnetic beam may be moving relative to the skin, the incident angle may be oblique, or the skin material properties may be inhomogeneous. Under these complex conditions, we no longer have a closed-form analytical solution for skin temperature, which needs to be solved numerically on a 3D grid and in general is not axially symmetric [18] [19]. Our volume computation methods are well capable of accommodating this complex situation. It requires only a discrete representation of skin temperature on a 3D grid; no symmetry is assumed. Here we use (26) to construct a non-symmetric problem for testing the accuracy of volume computation. We first find the activated volume of (26). The activated skin region of (26) is

$\begin{array}{c}D=\left\{\left(x,y,z\right)|P_d^{\left(0\right)}\exp \left(\frac{-\left(x^2+y^2\right)}{2}\right)U\left(z,t\right)\ge 1\right\}\\ =\left\{\left(x,y,z\right)|\left(x^2+y^2\right)\le 2\log \left(P_d^{\left(0\right)}U\left(z,t\right)\right),0\le z\le z_L\left(t\right)\right\}\end{array}$ (27)

where $z_L\left(t\right)\equiv U^{-1}\left(\frac{1}{P_d^{\left(0\right)}},t\right)$ and $U^{-1}\left(\cdot ,t\right)$ is the inverse of $z\mapsto U\left(z,t\right)$ at fixed $t$ . We express the activated skin volume in terms of a single integral.

$\text{volume}\left(D\right)=2\pi {\displaystyle {\int }_0^{z_L\left(t\right)}\text{log}\left(P_d^{\left(0\right)}U\left(z,t\right)\right)\text{d}z}$

We solve this integral accurately to machine precision using the composite Simpson’s rule, and use the result to assess the errors in numerical volumes (6) and (9). In the numerical test, we use $P_d^{\left(0\right)}=2.5$ and $t=1$ . To test the robustness of the volume computation, we shift the region in the x-direction by a z-dependent amount. This transformation destroys the axial symmetry while preserving the true volume, which is needed in error assessment. We set the test problem as follows.

$\{\begin{array}{l}\text{function}:f\left(x\right)=P_d^{\left(0\right)}\exp \left(\frac{-\left({\left(x+0.5z\right)}^2+y^2\right)}{2}\right)U\left(z,t\right)-1,x=\left(x,y,z\right)\\ \text{enclosed region}:D=\left\{\left(x,y,z\right)|f\left(x\right)\ge 0\right\}\\ \text{volume}\left(D\right)=1.29409605751414\left(\text{at}{P}_d^{\left(0\right)}=2.5,t=1\right)\end{array}$ (28)

Figure 10 shows a triangulation of the isosurface at the activation temperature threshold (left panel) and the errors in numerical volumes (right panel). The triangulation is based on a coarse grid with $\text{Δ}x=\frac{1}{8}$ for illustrative purposes. The activated skin region $D$ is enclosed by the isosurface and the skin surface ($z=0$ ). With the z-dependent shift in $x$ , neither function $f\left(x\right)$ in (28) nor the region $D$ enclosed by the isosurface is axially symmetric. The volume computation procedure, however, is not impacted at all by the transformation. It uses only a discretization of function $f\left(x\right)$ on a 3D numerical grid. The errors in numerical volumes behave as expected. The error in $V\left(\text{Δ}x\right)$ is proportional to ${\left(\text{Δ}x\right)}^2$ ; the error in $V_{\text{extrap}}\left(\text{Δ}x\right)$ is significantly below that in $V\left(\text{Δ}x\right)$ .

5. Summary

In this study, we investigated an accurate and reliable method for computing the activated skin volume in electromagnetic heating when skin is exposed to an incident electromagnetic beam of millimeter wavelength. The activated volume is enclosed by the isosurface of skin temperature at the level of thermal nociceptor activation threshold. Thus, the task of computing the activated volume is cast into a general problem of finding the volume enclosed by an isosurface of an underlying function $f\left(x\right)$ . Mathematically, the enclosed volume is expressed as a surface integral over the isosurface.

(a) (b)

Figure 10. Left: Triangulation of the activated skin region as an isosurface of $f\left(x\right)$ in (28). $f\left(x\right)$ has been transformed to destroy axial symmetry. Right: errors in numerical volumes $V\left(\text{Δ}x\right)$ and $V_{\text{extra}}\left(\text{Δ}x\right)$ vs mesh spacing $\text{Δ}x$ .

In many electromagnetic heating applications, the skin temperature is obtained numerically by solving a heat transfer model on a 3D grid. That is, the skin temperature is known only on a 3D grid. With only a discrete representation of function $f\left(x\right)$ , the corresponding isosurface is obtained approximately using linear interpolation along grid lines of the 3D grid. The resulting approximation of the isosurface is a collection of planar triangle elements. With a triangulation approximation of the isosurface, surface integration is carried out over each planar triangle element, which is analytically given by a triple scalar product. The enclosed volume is obtained by summing over all triple scalar products, each for one planar triangle element of the isosurface.

The process of approximating a surface using flat triangle elements is analogous to that of approximating a curve using line segments. Both are expected to produce second-order approximations. Given a discrete representation of $f\left(x\right)$ with mesh spacing $\text{Δ}x$ , a triangulation of the isosurface is constructed, and the enclosed volume is calculated. The dominant error in the numerical volume obtained with mesh spacing $\text{Δ}x$ is proportional to ${\left(\text{Δ}x\right)}^2$ . From the given fine representation of $f\left(x\right)$ , a coarse representation with mesh spacing $\left(2\text{Δ}x\right)$ is obtained by down-sampling, which does not involve any additional numerical simulations of a 3D heat transfer model. The dominant error in the numerical volume obtained with mesh spacing $\left(2\text{Δ}x\right)$ is proportional to ${\left(2\text{Δ}x\right)}^2$ . Extrapolation combines these two numerical volumes to eliminate the dominant error. The result of extrapolation is a more accurate approximation to the enclosed volume.

We evaluated the performance of the two volume computation methods: 1) triangulation only, 2) triangulation plus extrapolation. To accurately assess the errors in numerical volumes, we selected test problems that have closed-form analytical solutions, including a simplified case of skin electromagnetic heating. In all problems tested, we observed that 1) the error in the numerical volume directly from triangulation of isosurface is proportional to ${\left(\text{Δ}x\right)}^2$ where $\text{Δ}x$ is the mesh spacing of the 3D grid in the given discrete representation of function $f\left(x\right)$ ; and 2) the error in the numerical volume from triangulation plus extrapolation is much smaller than that in the triangulation-only volume, clearly demonstrating that extrapolation consistently improves the accuracy across all test cases.

Acknowledgement

The authors acknowledge the Joint Intermediate Force Capabilities Office of the U.S. Department of Defense and the Naval Postgraduate School for supporting this work.

Disclaimer

The views expressed in this document are those of the authors and do not reflect the official policy or position of the Department of Defense or the U.S. Government.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Romanenko, S., Begley, R., Harvey, A.R., Hool, L. and Wallace, V.P. (2017) The Interaction between Electromagnetic Fields at Megahertz, Gigahertz and Terahertz Frequencies with Cells, Tissues and Organisms: Risks and Potential. Journal of the Royal Society Interface, 14, Article ID: 20170585.[CrossRef] [PubMed]
[2]	Zhadobov, M., Chahat, N., Sauleau, R., Le Quement, C. and Le Drean, Y. (2011) Millimeter-Wave Interactions with the Human Body: State of Knowledge and Recent Advances. International Journal of Microwave and Wireless Technologies, 3, 237-247.[CrossRef]
[3]	Le Quément, C., Nicolaz, C.N., Habauzit, D., Zhadobov, M., Sauleau, R. and Le Dréan, Y. (2014) Impact of 60-GHz Millimeter Waves and Corresponding Heat Effect on Endoplasmic Reticulum Stress Sensor Gene Expression. Bioelectromagnetics, 35, 444-451.[CrossRef] [PubMed]
[4]	Mason, P.A., Walters, T.J., Nelson, M.T. and Nelson, D.A. (2000) Skin Heating Effects of Millimeter-Wave Irradiation-Thermal Modeling Results. IEEE Transactions on Microwave Theory and Techniques, 48, 2111-2120.[CrossRef]
[5]	Foster, K.R., Zhang, H. and Osepchuk, J.M. (2010) Thermal Response of Tissues to Millimeter Waves: Implications for Setting Exposure Guidelines. Health Physics, 99, 806-810.[CrossRef] [PubMed]
[6]	Cazares, S., Snyder, J.A., Belanich, J., Biddle, J., Buytendyk, A., Teng, S.H., et al. (2019) Active Denial Technology Computational Human Effects End-to-End Hypermodel (ADT CHEETEH). Human Factors and Mechanical Engineering for Defense and Safety, 3, Article No. 13. [Google Scholar] [CrossRef]
[7]	Topfer, F. and Oberhammer, J. (2015) Millimeter-Wave Tissue Diagnosis: The Most Promising Fields for Medical Applications. IEEE Microwave Magazine, 16, 97-113.[CrossRef]
[8]	Srivastava, A., Gupta, M.S. and Kaur, G. (2020) Energy Efficient Transmission Trends Towards Future Green Cognitive Radio Networks (5G): Progress, Taxonomy and Open Challenges. Journal of Network and Computer Applications, 168, Article ID: 102760.[CrossRef]
[9]	Sheen, D.M., McMakin, D.L. and Hall, T.E. (2007) Detection of Explosives by Millimeter-Wave Imaging. In: Yinon, J., Ed., Counterterrorist Detection Techniques of Explosives, Elsevier, 237-277.[CrossRef]
[10]	Liu, H. (2020) Robot Systems for Rail Transit Applications. Elsevier.
[11]	Wang, H., Burgei, W.A. and Zhou, H. (2020) A Concise Model and Analysis for Heat-Induced Withdrawal Reflex Caused by Millimeter Wave Radiation. American Journal of Operations Research, 10, 31-81.[CrossRef]
[12]	Wang, H., Burgei, W.A. and Zhou, H. (2020) Non-Dimensional Analysis of Thermal Effect on Skin Exposure to an Electromagnetic Beam. American Journal of Operations Research, 10, 147-162.[CrossRef]
[13]	Wang, H., Foley, S.E. and Zhou, H. (2025) Asymptotic Solution for Skin Heating by an Electromagnetic Beam at an Incident Angle. Electronics, 14, Article No. 3061.[CrossRef]
[14]	Skoog, D.A., Holler, F.J. and Crouch, S.R. (2017) Principal of Instrumental Analysis. 7th Edition, Sunder College Publisher.
[15]	Jaime-Yepez, U., Wang, H., Foley, S.E. and Zhou, H. (2024) Asymptotic Solution of Electromagnetic Heating of Skin Tissue with Lateral Heat Conduction. Journal of Engineering Mathematics, 147, Article No. 14.[CrossRef]
[16]	Walters, T.J., Blick, D.W., Johnson, L.R., Adair, E.R. and Foster, K.R. (2000) Heating and Pain Sensation Produced in Human Skin by Millimeter Waves: Comparison to a Simple Thermal Model. Health Physics, 78, 259-267.[CrossRef] [PubMed]
[17]	Wang, H., Burgei, W.A. and Zhou, H. (2020) Analytical Solution of One-Dimensional Pennes’ Bioheat Equation. Open Physics, 18, 1084-1092.[CrossRef]
[18]	Rodríguez-Pérez, A.S., Gilardi-Velázquez, H.E. and Velázquez-Pérez, S.E. (2025) Analysis and Simulation of Dynamic Heat Transfer and Thermal Distribution in Burns with Multilayer Models Using Finite Volumes. Dynamics, 5, Article No. 41.[CrossRef]
[19]	Joshi, A., Wang, F., Kang, Z., Yang, B. and Zhao, D. (2022) A Three-Dimensional Thermoregulatory Model for Predicting Human Thermophysiological Responses in Various Thermal Environments. Building and Environment, 207, Article ID: 108506.[CrossRef]