*Article* **Calculation of Effective Thermal Conductivity for Human Skin Using the Fractal Monte Carlo Method**

**Guillermo Rojas-Altamirano <sup>1</sup> , René O. Vargas 1,\* , Juan P. Escandón <sup>1</sup> , Rubén Mil-Martínez <sup>2</sup> and Alan Rojas-Montero <sup>1</sup>**


**Abstract:** In this work, an effective thermal conductivity (ETC) for living tissues, which directly affects the energy transport process, is determined. The fractal scaling and Monte Carlo methods are used to describe the tissue as a porous medium, and blood is considered a Newtonian and non-Newtonian fluid for comparative and analytical purposes. The effect of the principal variables—such as fractal dimensions *D<sup>T</sup>* and *D<sup>f</sup>* , porosity, and the power-law index, *n*—on the temperature profiles as a function of time and tissue depth, for one- and three-layer tissues, besides temperature distribution, are presented. ETC was improved by considering high tissue porosity, low tortuosity, and shearthinning fluids. In three-layer tissues with different porosities, perfusion with a non-Newtonian fluid contributes to the understanding of the heat transfer process in some parts of the human body.

**Keywords:** effective thermal conductivity; fractal scaling; Monte Carlo; porous media; non-Newtonian fluid; power-law model; bioheat equation; human body

## **1. Introduction**

The skin is the largest single organ of the body, enabling protection from the surrounding environment. It consists of several layers and plays an important role in thermoregulation, sensory, and host defense functions [1–3]. The skin is generally described by a three-layer tissue: epidermis, dermis, and hypodermis (also called subcutaneous) [4–6]. The thickness of these layers varies depending on the location of the skin. The epidermis is the outer layer (75–150 µm), this layer plays a barrier role between environment and organism [7,8]. The dermis is much thicker than the epidermis, in this, there are blood vessels, nerves, lymph vessels, and skin appendages. Dermis performances important functions in thermoregulation and supports the vascular network to supply the non-vascularized epidermis with nutrients. This layer is formed by an irregular network with wavy and unaligned collagen fiber bundles, allowing considerable deformations in all directions. The hypodermis is composed of loose fatty connective tissue. It is not part of the skin, but appears as a deep extension of the dermis, and depends on the age, sex, race, endocrine, and nutritional status of the individual [7–10]. The thermoregulation function of skin is realized mainly by modifying the blood flow, which is located in a microcirculatory bed, composed of arterioles, arterial and vein capillaries, and venules (blood perfusion). Blood perfusion has great effect on the heat transfer process in living tissues [8,11,12]. Heat transfer in human tissues takes place through different mechanisms, such as heat conduction, blood perfusion, metabolic heat generation, and external interactions [13,14]. One of the earliest models of heat transfer in biological tissues was developed by Pennes in 1948 [15], who proposed a model to describe the effects of metabolism and blood perfusion on the energy balance within tissue, this model is based on the classical Fourier's law [15]. Wulff [16]

**Citation:** Rojas-Altamirano, G.; Vargas, R.O.; Escandón, J.P.; Mil-Martínez, R.; Rojas-Montero, A. Calculation of Effective Thermal Conductivity for Human Skin Using the Fractal Monte Carlo Method. *Micromachines* **2022**, *13*, 424. https:// doi.org/10.3390/mi13030424

Academic Editor: Lanju Mei and Shizhi Qian

Received: 4 February 2022 Accepted: 6 March 2022 Published: 10 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

questioned the assumptions of the Pennes model and provided an alternative analysis. He assumed that heat transfer between blood flow and tissue should be modeled proportionally to the temperature difference between these two media and not between the two temperatures of the blood flow. Klinger [17] consider the convective heat transfer caused by the blood flow inside the tissue, since this term was neglected by Pennes. Chen and Holmes [18] assumed that the total tissue control volume is composed of the solid-tissue subvolume and blood subvolume. They determined an effective thermal conductivity using the tissue porosity and the local mean tissue temperature, together with a simplified volume-averaging technique for the solid and tissue spaces. Weinbaum et al. [19,20] determined an effective thermal conductivity (ETC) based on the hypothesis that small arteries and veins are parallel and the flow direction is countercurrent, which is a function of the blood flow rate and vascular geometry. The model included a perfusion bleed-off term that apparently resembles the Pennes perfusion term. Weinbaum and Jiji [11] derived a simplified equation to study the influence of the blood flow on the tissue temperature distribution defining an ETC.

In biological tissues, many body parts reveal anisotropy in heat transport that can not be explained by the Fourier's law [6,8]. This leads to formulation of thermal wave bioheat model based on two main approaches: the Maxwell–Cattaneo approach with heat flux time lag (also known as the single-phase approach), and the double-phase-lag (DPL) approach with relaxations in both the heat flux and temperature gradient propagation [21–23]. Various researchers have contributed in this area with analytical and experimental work. Hobiny and Abbas [24] presented an analytical solution of the hyperbolic bioheat equation under intense moving heat source. Alzahrani and Abbas [25] also presented an analytical approach, experimental temperature data, and a time sequential concept to obtain the thermal damage and temperature in a living tissue due to laser irradiation. Hobiny and Abbas [26] provided a method to determine numerical solutions for thermal damage of cylindrical living tissues using hyperbolic bioheat model. Hobiny et al. [27] proposed a new interpretation to study thermal damage in a skin tissue caused by laser irradiation, using the fractional order bioheat model. Hobiny et al. [28] presented an analytical method and experimental verification, to estimate thermal damage and temperature due to laser irradiation, using skin surface measurement data. Kumari and Singh [29] generated a spacefractional mathematical model of bioheat transfer to graphically analyze thermal behavior within living tissue, using a three-phase-lag constitutive relation. Li et al. [30] developed a generalized model of bioheat transfer to explore heat transport properties involving different thermal phase lagging effect. Important reviews and articles that the reader can consult additionally in this context are the following: classical mathematical models of bioheat [13]; developments in modeling heat transfer in blood perfused tissues [9,31]; bioheat models based on the porous media theory [32]; general heat transfer review [33,34]; concepts, derivation, and experimental versus porous media modeling [35]; and modeling and scaling of the bioheat equation [3].

Alternatively to continuum models, concepts that consider the tissue matrix, arteries, veins, and capillary vessels in a porous medium with specific porosity variations, ETC, and heat dispersion by blood flow have been developed [3,12,13,32,36]. Porous medium is defined as a material volume consisting of solid matrix with an interconnected pores. It is characterized by porosity, ratio of the pore space to the total volume of the medium, permeability, and tortuosity [32]. Khanafer and Vafai [36] remarked that the most appropriate treatment for heat transfer in biological tissues is the porous media theory because of fewer assumptions as compared with other models. Roetzel and Xuan [37] introduced a two-equation bioheat model in which the biological system is a porous media. It is divided into two different regions, the vascular and extravascular, without considering local thermal equilibrium between the two phases, introducing an equivalent ETC in the energy equations of blood and tissue [37]. Nakayama and Kuwahara [38] developed a model that consists of two energy equations based on the volume average theory (VAT), these equations are correct for all cases of thermal non-equilibrium.

The ETC is one of the most important thermo-physical properties for quantifying conductive heat transfer of porous media with gas, liquid, and solid phases [39]. The prediction of the ETC of porous media is essential to many engineering applications, such as thermally enhanced oil recovery, geothermal energy, and chemical and biological engineering [40]. With the development of computer technology, many numerical methods, such as Monte Carlo [41] and lattice Boltzmann [42], have been proposed to study the conductive heat transfer and evaluate the ETC of porous media. In addition to conventional methods based on Euclidean geometry, fractal geometry has been shown with evident advantages for addressing the complexity and multiple scales of porous media [43]. Hence, the fractal geometry has been successfully applied to characterize structures of transport processes in porous media [44,45]. Kou et al. proposed a fractal model for ETC of porous media based on fractal scaling law for water and gas phases in the pores [46]. Extensive studies have shown that most natural porous media and some synthetic porous media possess self-similar fractal scaling laws over multiple scales [45]. Therefore, two kinds of fractal models based on pore and solid phases have been proposed [47]. Fractal scaling laws can be applied to characterize the geometrical and morphological structures for pore and solid phases in porous media, respectively.

The fractal Monte Carlo method has been applied in different areas to model a porous media. Yu et al. [48] performed Monte Carlo simulations to predict the permeability of fractal porous media, their results were verified by comparison with the analytical solution for the permeability of bi-dispersed porous media. Zou et al. [49] used the Monte Carlo simulation technique to model the surface topography in a scale-invariant manner with the fractal nature of rough surfaces. Yu [44] presented a review article summarizing the theories, methods, mathematical models, achievements, and open questions in the area of flow in fractal porous media by applying the theory and technique of fractal geometry. Feng et al. [50] combined the Monte Carlo technique with fractal geometry theory to predict the thermal conductivity of nanofluids. Xu et al. [51] performed Monte Carlo simulations of radial seepage flow in the fractured porous medium, where the fractal probability model was applied to characterize the fracture size distribution. Vadapalli et al. [52] proposed a permeability estimation method for a sandstone reservoir, which considers the fractal behavior of pore size distribution and tortuosity of capillary pathways using Monte Carlo simulations. Xu et al. [53] used fractal Monte Carlo simulations to predict the effective thermal conductivity of porous media. Xiao et al. [54] employed the fractal Monte Carlo to simulate the Kozeny–Carman constant of fibrous porous media with the micropore size characterized by the fractal scaling law. Yang et al. [55] performed Monte Carlo simulations based on the fractal probability law to understand gas flow mechanisms and predict the apparent gas permeability of shale reservoirs.

In this work, the calculation of ETC for human skin using the fractal scaling and Monte Carlo methods is presented. This ETC involves a bundle of tortuous capillaries whose size distribution follow fractal scaling laws. The power-law model was chosen because of its simplicity in describing different non-Newtonian fluids by modifying a single parameter, in the case of blood as a shear-thinning fluid. The heat transfer process in the perfused tissue is analyzed, and a heat source is applied on the tissue surface for a period of time without reaching the degradation temperature. The tissue is considered as a uniform porous medium of one and three layers, which can be assigned different porosity and conductivity. The temperature profiles, as well as their distribution when modifying the main variables of the model, are presented. To the authors' knowledge, there are no studies that determine an ETC for heat transfer in biological tissues, considering the techniques mentioned above and especially for non-Newtonian fluids.

#### **2. Heat Transfer in Human Skin**

Figure 1 presents the human skin structure, considering three layers; epidermis, dermis and hypodermis. These layers differ in having their own physical properties such as density, specific heat, thermal conductivity and porosity.

**Figure 1.** Human skin structure.

The first equation that described heat transfer in human tissue and included the effects of blood flow on tissue temperature on a continuum basis was presented by Pennes [15]. He presented the heat transfer analysis in the human forearm, considering the metabolic heat rate in the tissue and the perfusion heat source term. This last term has been the focus of attention since its inception and many researches have found alternative representations of the effect of blood perfusion on tissue heat transfer. The Pennes equation is given by:

$$
\rho\_l c\_l \frac{\partial T\_l}{\partial t} = \nabla \cdot (k\_l \nabla T\_l) + \rho\_b c\_b \omega\_b (T\_a - T\_l) + q\_{m\nu} \tag{1}
$$

where *ρ*, *c*, *T*, *t*, *k*, *qm*, and *ω<sup>b</sup>* are density, specific heat, temperature, time, thermal conductivity, metabolic heat production, and blood perfusion rate per unit volume of tissue, respectively. The subscripts *t*, *b*, and *a* refer to tissue, blood, and artery, respectively [15] .

However, some inconsistencies in the Pennes model include the following: the thermal equilibrium take place in arteries and veins (not in the capillaries, as it assumes); it does not take into account any vascular architecture; and the most critical assumption is on the blood perfusion term, which is not a global term—it is local along the capillary and depends on direction.

#### **3. Mathematical Modeling**

Hyperthermia treatment consists of applying heat in a specific area of the human body. In this work, an external heat source is applied to an area of the forearm, as shown in Figure 2a. In addition, in Figure 2b, it is described that *H* is the total thickness of the tissue and the thicknesses of each layer are *H<sup>E</sup>* = 0.04*H*, *H<sup>D</sup>* = 0.48*H*, and *H<sup>H</sup>* = 0.48*H*. Initially, the tissue is at a constant temperature, *T<sup>c</sup>* ∼ 37 ◦C, subsequently, a heat source is applied to an area of the tissue. The area around the heat application area is open to the surroundings, and generally is a temperature lower than the human body. For this work, the surrounding temperature is considered to be *T*<sup>∞</sup> ∼ 25 ◦C. Furthermore, Figure 2b shows a mathematical representation of the hyperthermia treatment, note that the deepest internal temperature is maintained at the body temperature. According to Weinbaum and Jiji [11], the simplified two-dimensional governing equation of heat transport in biological tissue is given by:

$$(\rho \mathbb{C}\_p)\_t \frac{\partial T}{\partial t} = \frac{\partial}{\partial \mathbf{x}} \left( k\_{eff} \frac{\partial T}{\partial \mathbf{x}} \right) + \frac{\partial}{\partial y} \left( k\_{eff} \frac{\partial T}{\partial y} \right) + q\_{m\nu} \tag{2}$$

where *ke f f* is the effective thermal conductivity (ETC) and *q<sup>m</sup>* is the metabolic heat.

According to physical model showed in Figure 2, Equation (2) is subject to the following initial and boundary conditions:

$$T(\mathfrak{x}, y, 0) = T\_{\mathfrak{c}\prime} \tag{3}$$

$$\frac{\partial}{\partial \mathbf{x}} T(0, y, t) = 0,\tag{4}$$

$$\frac{\partial}{\partial \mathbf{x}} T(\mathcal{W}, y, t) = 0,\tag{5}$$

$$T(\mathbf{x}, \mathbf{0}, t) = T\_{\mathbf{c}} \tag{6}$$

$$-k\_l \frac{\partial}{\partial y} T(\mathbf{x}, H, t) = h[T\_\infty - T(\mathbf{x}, H, t)], \quad \frac{2}{5}W \ge \mathbf{x} \ge \frac{3}{5}W,\tag{7}$$

$$-k\_{\rm f} \frac{\partial}{\partial y} T(\mathbf{x}, H, t) = f\_{\prime} \quad \frac{2}{5} \mathcal{W} < \mathbf{x} > \frac{3}{5} \mathcal{W}\_{\prime} \tag{8}$$

$$f = \begin{cases} \begin{array}{c} q\_{\rm app} = 200 \left[ \frac{\rm W}{\rm m^2} \right]\_{\prime} \qquad t \le t\_{\rm app} \\\ h \big[ T\_{\infty} - T(x\_{\prime} H\_{\prime} t) \big]\_{\prime} \quad t > t\_{\rm app} \end{array} \end{cases}$$

where *W* and *H* are the width and the height of the domain, respectively. The *h*, *q*app, and *t*app are the heat transfer coefficient, applied external heat, and application time, respectively.

**Figure 2.** (**a**) Hyperthermia treatment; (**b**) human skin three-layer model, and the corresponding boundary conditions.

ETC is an important parameter in Equation (2). Weinbaum and Jiji [11] proposed a vascular function *V*(*y*) that can be constructed knowing the vascular data, which is the distribution of the arteries, veins, and capillaries. According to Weinbaum and Jiji [11], the vascular function increases with tissue depth [14]. In this work, the *representative elementary volume* (REV), is defined and the fractal scaling method is used to depict the vascular tissue structure.

## *3.1. Fractal Scaling Method*

We consider a cubic REV as shown in Figure 3, with defined length side *L*0. All REV capillaries extend throughout the volume from one side to the other as is showed in Figure 3.

**Figure 3.** Representative elementary volume of human skin.

The fractal scaling method establishes the relationship between the number and pore size in the porous medium. The fundamental fractal scaling law is applied to REV crosssection, as follows [56]:

$$N(>\lambda) = \left(\frac{\lambda\_{\text{max}}}{\lambda}\right)^{D\_f} \tag{9}$$

where *N*, *D<sup>f</sup>* , *λ*, and *λ*max are the number of capillaries, the fractal dimension, the equivalent diameter, and the maximum equivalent diameter of the capillaries in the REV cross-section, respectively. The number of capillaries with equivalent diameter between *λ* + *dλ* in the REV cross-section is:

$$-d\mathcal{N}(\lambda) = D\_f \lambda\_{\text{max}}^{D\_f} \lambda^{-D\_f - 1} d\lambda. \tag{10}$$

The total number of pores in the range from *λ*min to *λ*max is obtained using Equation (9), as follows:

$$N\_T(\leq \lambda\_{\min}) = \left(\frac{\lambda\_{\max}}{\lambda\_{\max}}\right)^{D\_f} \tag{11}$$

where *N<sup>T</sup>* is the total pores. Dividing Equation (10) by Equation (11), we obtain:

$$-\frac{dN}{N\_T} = D\_f \lambda\_{\text{min}}^{D\_f} \lambda^{-\binom{D\_f+1}{2}} d\lambda = f(\lambda) d\lambda,\tag{12}$$

where *f*(*λ*) is the probability density function, which satisfies that *f*(*λ*) ≥ 0. According to probability theory, *f*(*λ*) must satisfy the following normalization relation or total cumulative probability:

$$-\int\_{\lambda\_{\rm min}}^{\lambda\_{\rm max}} \frac{dN}{N\_T} = \int\_{\lambda\_{\rm min}}^{\lambda\_{\rm max}} f(\lambda) d\lambda = 1 - \left(\frac{\lambda\_{\rm min}}{\lambda\_{\rm max}}\right)^{D\_f} \equiv 1. \tag{13}$$

The integration result of Equation (13) shows that it holds if—and only if—the following holds:

$$\left(\frac{\lambda\_{\text{min}}}{\lambda\_{\text{max}}}\right)^{D\_f} \cong 0.\tag{14}$$

The above equation implies that *λ*min *λ*max, it must be satisfied for fractal analysis of a porous media. According to Yu et al. [48], Equation (14) can be considered as a criterion whether a porous medium can be characterized by fractal theory and technique. If Equation (13) is expressed as:

$$R(\lambda) = \int\_{\lambda\_{\rm min}}^{\lambda} f(\lambda) d\lambda = 1 - \left(\frac{\lambda\_{\rm min}}{\lambda}\right)^{D\_f} \tag{15}$$

then the pore diameter *λ* can be found as:

$$
\lambda = \frac{\lambda\_{\rm min}}{\left(1 - R\right)^{1/D\_f}} = \left(\frac{\lambda\_{\rm min}}{\lambda\_{\rm max}}\right) \frac{\lambda\_{\rm max}}{\left(1 - R\right)^{1/D\_f}}.\tag{16}
$$

On the other hand, the fractal path of the tortuous capillary can be described as follows [57]:

$$L(\lambda) = L\_0^{D\_T} \lambda^{1 - D\_T} \,. \tag{17}$$

where *L*(*λ*), *L*0, and *D<sup>T</sup>* are capillary tortuous length, characteristic length of a straight capillary (REV side length), and fractal dimension describing the capillary tortuous length, respectively. An important parameter involved in the ETC is the total pore area, *Ac*, determined by Wu and Yu [58]:

$$A\_{\mathcal{E}} = -\int\_{\lambda\_{\rm min}}^{\lambda\_{\rm max}} \frac{\pi \lambda^2}{4} dN(\lambda) = \frac{\pi D\_f \lambda\_{\rm max}^{D\_f} \left(\lambda\_{\rm max}^{2-D\_f} - \lambda\_{\rm min}^{2-D\_f}\right)}{4\left(2 - D\_f\right)},\tag{18}$$

the total cross area, *AT*, is:

$$A\_T = L\_0^2 = \frac{A\_c}{\Phi} \,\prime \tag{19}$$

where *φ* is the porosity.

#### *3.2. The Fractal ETC of the REV*

In order to obtain the ETC and according to Fourier's law, the total heat flux in the REV is given by:

$$Q\_T = k\_{eff} A\_T \frac{\Delta T}{L\_0} \,\tag{20}$$

where *ke f f* is the ETC including the tissue conductivity and convective blood flow in the capillaries, and ∆*T* is the difference temperature between two faces in the REV. The heat flux in a single tortuous capillary of the REV is:

$$Q\_c = k\_b \frac{\pi \lambda^2}{4} \frac{\Delta T}{L(\lambda)} = k\_b \frac{\pi \lambda^2}{4} \frac{\Delta T}{L\_0^{D\_T} \lambda^{1 - D\_T}} \, \tag{21}$$

where *k<sup>b</sup>* is the blood thermal conductivity. The heat flux corresponding only to the tissue is expressed by:

$$Q\_t = k\_t A\_l \frac{\Delta T}{L\_0} = k\_l (1 - \phi) A\_T \frac{\Delta T}{L\_0} \,\tag{22}$$

where the subscript *t* refers to the tissue. According to superposition theorem, the total heat flux *Q<sup>T</sup>* is the addition of the fluxes as follows:

$$Q\_T = Q\_l + Q\_{\mathcal{L}}.\tag{23}$$

Substituting Equations (21) and (22) into Equation (23) the following equation is obtained:

$$k\_{eff} = k\_l(1 - \phi) + k\_b \frac{\phi \left(2 - D\_f\right) \lambda^{D\_T + 1}}{L\_0^{D\_T - 1} D\_f \lambda\_{\text{max}}^{D\_f} \left(\lambda\_{\text{max}}^{2 - D\_f} - \lambda\_{\text{min}}^{2 - D\_f}\right)},\tag{24}$$

where *L*<sup>0</sup> is determined from Equation (19). According to the work presented by Weinbaum and Jiji [11], from an energy conservation balance in a countercurrent artery and mean value theory, derived the effective enhancement of the tissue conductivity. Therefore, by comparing the Weinbaum and Jiji ETC to Equation (24), it can be seen that the two equations have a similar form, proposing the following the expression:

$$k\_{eff} = k\_t (1 + \Gamma(\mathbf{r})\psi(\mathbf{r})),\tag{25}$$

where Γ is a parameter that depends on conduction or convective heat transport [59], *ψ* is related to skin structure (also called dimensionless vascular geometry function) [14], and **r** is the position vector. This work focuses on the convective transport of blood. Hence, the dimensionless ETC for this work is given by:

$$k\_{eff}^{\*} = \frac{k\_{eff}}{k\_l} = \left[ (1 - \phi) + \mathrm{Pe}^2 \left( \frac{k\_b}{k\_l} \right)^2 \frac{\phi \left( 2 - D\_f \right) \lambda^{D\_T + 1}}{L\_0^{D\_T - 1} D\_f \lambda\_{\mathrm{max}}^{D\_f} \left( \lambda\_{\mathrm{max}}^{2 - D\_f} - \lambda\_{\mathrm{min}}^{2 - D\_f} \right)} \right], \tag{26}$$

where *k* ∗ *e f f* is the dimensionless ETC and *Pe* is the Peclet number defined as Pe = PrRe = *ρbCbλu*/*k<sup>b</sup>* ; where *ρ<sup>b</sup>* , *C<sup>b</sup>* , and *u* are the density, specific heat, and average blood velocity in the capillary, respectively.

## *3.3. Fractal Dimensions, D<sup>f</sup> and D<sup>T</sup>*

Equation (26) is a function of porosity, fractal dimensions, maximum and minimum diameters, conductivity ratio, and Peclet number, the latter being a function of physical properties, velocity, and diameter.

There is a relationship between porosity and fractal dimension, *D<sup>f</sup>* , according to Yu and Li [56], is given by:

$$
\phi = \left(\frac{\lambda\_{\rm min}}{\lambda\_{\rm max}}\right)^{D\_E - D\_f} \tag{27}
$$

where *D<sup>E</sup>* is the Euclidean dimension (*D<sup>E</sup>* = 2 and 3, for two- and three-dimensional space, respectively). Another important aspect of Equation (27), from experimentation, it has been found that the ratio of minimum and maximum diameter, *λ*min/*λ*max, in several natural porous media, is the order of <sup>10</sup>−<sup>2</sup> <sup>∼</sup> <sup>10</sup>−<sup>4</sup> , [60]. Feng et al. derived a generalized model covering a wide range of porosities, for the effective thermal conductivity, based on the fact that statistical self-similarity exists in porous media [61].

In this work, *D<sup>T</sup>* is established manually, taking into account whether *D<sup>T</sup>* > 1 the tortuosity is present and *D<sup>T</sup>* = 1 are straight capillaries.

#### *3.4. Non-Newtonian Fractal Velocity*

The Peclet number in Equation (26) is defined by:

$$\text{Pe} = \text{PrRe} = \frac{\rho\_b \mathcal{C}\_b}{k\_b} \lambda u\_\prime \tag{28}$$

where the velocity *u* is a function of the microcapillary diameter, is obtained from the non-Newtonian fluid flow in a single microcapillary, as presented by Zhang [62], as follows:

$$q\_c = \left[\frac{dp}{dL\_0} \frac{L\_0^{1-D\_T} 2^{D\_T-1}}{\mu\_b (D\_T+1) D\_T}\right]^{\frac{1}{n}} \frac{n\pi}{D\_T+3n} \left(\frac{\lambda}{2}\right)^{\frac{D\_T}{n}+3} \tag{29}$$

where *qc*, *p*, and *µ<sup>b</sup>* are the flow rate in a single microcapillary, pressure, and dynamic blood viscosity, respectively. By considering *q<sup>c</sup>* = *uAsc*, where *Asc* is the area of a single capillary and the velocity can be determined as:

$$u = \left[\frac{dp}{dL\_0} \frac{L\_0^{1-D\_T} 2^{D\_T-1}}{\mu (D\_T+1) D\_T}\right]^{\frac{1}{n}} \frac{n}{D\_T+3n} \left(\frac{\lambda}{2}\right)^{\frac{D\_T}{n}+1} \,\tag{30}$$

where *n* is the power-law index of the constitutive equation, when *n* = 1 the Newtonian velocity is recovered. For *n* < 1 and *n* > 1 describes shear-thinning and shear-thickening fluid behavior, respectively. Figure 4 presents the average velocity and Peclet number as a function of the capillary diameter for different values of the power-law index, *n*, considering two values of porosity, *φ* = 0.1 and *φ* = 0.5. By increasing the pore diameters, the flow through the tissue is greater for all cases. This effect is magnified when the power-law index, *n*, decreases, which indicates a lower resistance to flow, because the viscosity decreases, which is characteristic of shear-thinning fluids. The opposite case can be seen in this figure for shear-thickening fluids (*n* > 1). The number of Peclet is proportional to the velocity, therefore they have the same tendency. On the other hand, by increasing the porosity, both *u* and *Pe* have a slight increase because the pores are very small. Figure 5 shows the average velocity and Peclet number as a function of the fractal tortuosity, *DT*, for different values of porosity, *φ*, considering two values of power-law index, *n* = 1 and *n* = 0.6. For this case, the maximum pore diameter was taken into account. Keeping constant porosity and increasing tortuosity, *DT*, the velocity and Peclet number both decrease. This is correct for complex vascular architectures where the flow experiences higher resistance as well as being very small. For straight or slightly tortuous capillaries—where the highest velocities and Peclet numbers are found—there is a large increase as the power-law index decrease, as shown in Figure 5.

**Figure 4.** (**a**) Average velocity and (**b**) Peclet number as a function of capillary diameter for different values of the power-law index *n*. In a single microcapillary—continuous lines (*φ* = 0.1), and dashed lines (*φ* = 0.5).

**Figure 5.** (**a**) Average velocity and (**b**) Peclet number as a function of fractal tortuosity *D<sup>T</sup>* for different porosity values. For both in a single microcapillary, continuous lines (Newtonian fluid), and dashed lines (non-Newtonian fluid).

#### *3.5. Monte Carlo Method*

The Monte Carlo method is used to establish the pore size and the porous medium distribution, as is shown in Figure 6. Figure 6a,b present one- and three-layer tissue with different porosities, with *φ* = 0.1 (one layer) and *φ* = 0, *φ* = 0.5, and *φ* = 0.05 for epidermis, dermis, and hypodermis, respectively. Yu et al. [48] was the first to propose the fractal Monte Carlo methodology to simulate the transports in fractal porous media. In Figure 6c for one-layer tissue and Figure 6d for three-layer tissue show the Monte Carlo simulations, which are performed in the range of *λ*min–*λ*max. Figure 6d shows the pore sizes variation for dermis and hypodermis layers, since the porosity of the epidermal layer is not considered. The pore size variation is determined using the Monte Carlo method in Equation (16), as follows:

$$
\lambda\_i = \frac{\lambda\_{\rm min}}{\left(1 - R\_i\right)^{1/D\_f}} = \left(\frac{\lambda\_{\rm min}}{\lambda\_{\rm max}}\right) \frac{\lambda\_{\rm max}}{\left(1 - R\_i\right)^{1/D\_f}}\,\tag{31}
$$

where *λ<sup>i</sup>* is the variation of the pore diameter and *R<sup>i</sup>* are the random numbers between 0 − 1. On the other hand, the random numbers also help us to construct the porous medium presented in Figure 6. Equation (31) is derived from Equation (9), which implies that there is only one largest pore in the REV cross-section [44]. This is consistent to the pore size distribution shown in Figure 6.

#### *3.6. Dimensionless Governing Equation*

The dimensionless variables are:

$$
\mathfrak{x} = \frac{\mathfrak{x}}{H}; \qquad \mathfrak{y} = \frac{\mathfrak{y}}{H}; \qquad \mathfrak{x} = t\frac{\mathfrak{a}\_{\mathfrak{f}}}{H^2}; \qquad \theta = \frac{T - T\_{\infty}}{T\_c - T\_{\infty}}.\tag{32}
$$

where *α<sup>t</sup>* is the tissue thermal diffusivity. Substituting dimensionless variables of Equation (32) into Equation (2), the dimensionless governing equation is:

$$\frac{\partial\theta}{\partial\pi} = \frac{\partial}{\partial\mathfrak{x}} \left( k\_{\varepsilon f f}^{\*} \frac{\partial\theta}{\partial\mathfrak{x}} \right) + \frac{\partial}{\partial\mathfrak{y}} \left( k\_{\varepsilon f f}^{\*} \frac{\partial\theta}{\partial\mathfrak{y}} \right) + \Phi\_{m\nu} \tag{33}$$

where *k* ∗ *e f f* = *ke f f* /*k<sup>t</sup>* is defined in Equation (26) and <sup>Φ</sup>*<sup>m</sup>* <sup>=</sup> *<sup>q</sup>mH*2/*kt*(*T<sup>c</sup>* <sup>−</sup> *<sup>T</sup>*∞) is the dimensionless metabolic heat generation. Respectively, the dimensionless initial and boundary conditions are given by:

$$\theta(\mathfrak{x}, \mathfrak{y}, \mathbf{0}) = 1,\tag{34}$$

$$\frac{\partial}{\partial \mathfrak{X}} \theta(0, \mathfrak{J}, \mathfrak{r}) = 0,\tag{35}$$

$$\frac{\partial}{\partial \mathfrak{X}} \theta(\mathfrak{Z}, \mathfrak{Y}, \mathfrak{r}) = 0,\tag{36}$$

$$
\theta(\vec{x},0,\tau) = 1,\tag{37}
$$

$$\frac{\partial}{\partial \bar{y}} \theta(\bar{\mathbf{x}}, 1, \tau) = Bi[\theta\_{\infty} - \theta(\bar{\mathbf{x}}, 1, \tau)], \quad \frac{2}{5} \ge \bar{\mathbf{x}} \ge \frac{3}{5},\tag{38}$$

$$\frac{\partial}{\partial \overline{y}} \theta(\mathfrak{x}, 1, \mathfrak{r}) = f^\*, \quad \frac{2}{5} < \mathfrak{x} > \frac{3}{5} \,\, \text{s} \tag{39}$$

$$f^\* = \begin{cases} \frac{H}{k\_l \Delta T} q\_{\rm app} & \tau \le \tau\_{\rm app}, \\\ Bi[\theta\_\infty - \theta(\mathfrak{x}, 1, \tau)], & \tau > \tau\_{\rm app}. \end{cases}$$

where *Bi* is the Biot number defined as *Bi* = *hH*/*k<sup>t</sup>* , *θ*<sup>∞</sup> is the dimensionless temperature of the environment, and *f* ∗ takes the boundary condition according to the dimensionless time of application *τ*app.

#### **4. Results**

The numerical results presented in this section were generated from a numerical code developed in the Fortran programming language. Equation (33) subject to boundary conditions Equations (34)–(39) was solved using an explicit finite difference method.

The values used for the tissue physical properties are: *ρ<sup>t</sup>* = 1200 kg/m<sup>3</sup> , *Cp*,*<sup>t</sup>* = 3600 J/kg· ◦C, *k<sup>t</sup>* = 0.293 W/m· ◦C. When considering three-layer tissue, their corresponding thermal conductivities are: *k<sup>e</sup>* = 0.25 W/m· ◦C, *k<sup>d</sup>* = 0.45 W/m· ◦C, and *k<sup>h</sup>* = 0.2 W/m· ◦C, for epidermis, dermis, and hypodermis, respectively. The blood physical properties are: *ρ<sup>b</sup>* = 1052 kg/m<sup>3</sup> , *Cp*,*<sup>b</sup>* = 3800 J/kg· ◦C, and *k<sup>b</sup>* = 0.582 W/m· ◦C [2,63]. The metabolic heat is *q<sup>m</sup>* = 368.1 W/m<sup>3</sup> [63]. The heat transfer coefficient due to convection and radiation in surroundings is: *h* = 5 W/m<sup>2</sup> · ◦C. The environmental temperature was chosen as *T*<sup>∞</sup> = 25 ◦C. For the three-layer model, the main changes in porosity are considered to be in the dermis, due to the greater interaction between the tissue and the vascular network [8]. This work does not take into account thermal degradation (thermal damage) in tissue, which occurs at 44 ◦C and higher [64].

#### *4.1. ETC Analysis*

According to Equation (26), ETC depends on the conductivity ratio, porosity, *φ*, fractal coefficients, *D<sup>f</sup>* , *DT*, and Peclet number. In this work, both conductivities of tissue and blood are constant; therefore, the conductivity ratio is also constant. Once the porosity value is assigned manually, *D<sup>f</sup>* can be determined using Equation (27). The Peclet number depends on the physical properties of the blood, capillary diameter, and velocity, see Equation (28). An important contribution of this work is that it allows the analysis of blood as a Newtonian and non-Newtonian fluid. This is possible through the power-law model, which is immersed in the velocity calculation Equation (30). This equation also depends on the fractal coefficients, in addition to the viscosity and the power-law index, *n*. The latter makes it possible to consider a Newtonian fluid (*n* = 1) and non-Newtonian fluids such as shear-thinning fluid (*n* < 1) and shear-thickening *n* < 1. The blood is a shear-thinning fluid described by a power-law index value of *n* = 0.6, according to Johnston et al. [65].

Figure 7a presents the ETC as a function of the fractal coefficient *D<sup>T</sup>* for different porosity values, in addition for two values of the power-law index *n* = 1 and *n* = 0.6. *D<sup>T</sup>* describes the capillary tortuous length, the influence of *D<sup>T</sup>* on the ETC is in the range

of 1–1.3, which is in accordance with that reported by Yu and Cheng [57]. By increasing the porosity, the ETC also increases, due to higher blood flow through the tissue having a higher thermal conductivity. An important increase is found in the ETC when considering blood as a non-Newtonian fluid, since it increases with porosity as mentioned above, but especially for a shear-thinning fluid such as blood (*n* = 0.6).

Figure 7b shows the fractal dimension as a function of porosity. When porosity tends to unity, the fractal dimension tends to two. That for a surface indicates that it is totally occupied by pores, which corresponds to a fractal dimension of two, it is consistent according to Yu et al. [66].

**Figure 7.** (**a**) Effective thermal conductivity as a function of fractal dimension *DT*. Continuous lines (Newtonian fluid) and dashed lines (non-Newtonian fluid). (**b**) Fractal dimension *D<sup>f</sup>* as a function of porosity, for different ratios *λ*min/*λ*max.

Figure 8 exhibits the ETC as a function of fractal dimension *D<sup>f</sup>* , and the porosity *φ*; for two pore diameter values, *<sup>λ</sup>* <sup>=</sup> <sup>4</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> and *<sup>λ</sup>*max <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> , and by considering straight capillaries *D<sup>T</sup>* = 1 for different values of the power-law index. *k* ∗ *e f f* has an increase by augmenting the pore diameter, but the main improvement is for shear-thinning fluids, due to the lower resistance that the fluid experiences in large pores, increasing the velocity and generating a higher energy dissipation.

**Figure 8.** Effective thermal conductivity as a function of (**a**) fractal dimension *D<sup>f</sup>* , and (**b**) porosity *φ*. Continuous lines (Newtonian fluid) and dashed lines (non-Newtonian fluid).

#### *4.2. Code Validation*

To validate the ETC, the solution of the Weinbaum and Jiji [11] (WJ) is compared with the present work using a one-layer tissue with porosity *φ* = 0.05, straight capillaries *D<sup>T</sup>* = 1, *kb*/*k<sup>t</sup>* = 1.9863, *D<sup>f</sup>* is determined by Equation (27) and Newtonian fluid *n* = 1. The Peclet number depends on the pore diameter *λ* and the velocity *u*, WJ determined *Pe* = 20, which is in accordance with this work of *Pe* = 16.3 for maximum pore size. Figure 9 shows the comparison of the temperature profile as a function on time and a function of deep-tissue layer (at *τ* = 0.09) between WJ model and the present work. The ETC obtained by Weinbaum and Jiji [11] is as follows:

$$\frac{k\_{\ell f, Wf}}{k\_t} = \left[1 + P e\_0^2 V(\vec{y})\right]\_{\prime} \tag{40}$$

where *Pe*<sup>0</sup> = 2*ρ<sup>b</sup> cba*0*u*0/*k<sup>b</sup>* = 20, *V*(*y*¯) = *A* + *By*¯ + *Cy*¯ <sup>2</sup> with *<sup>A</sup>* <sup>=</sup> 6.32 <sup>×</sup> <sup>10</sup>−<sup>5</sup> , *B* = <sup>−</sup>15.9 <sup>×</sup> <sup>10</sup>−<sup>5</sup> and *<sup>C</sup>* <sup>=</sup> <sup>10</sup> <sup>×</sup> <sup>10</sup>−<sup>5</sup> . The metabolic heat Φ*<sup>m</sup>* = *qmL* <sup>2</sup>/*kt*(*T<sup>c</sup>* <sup>−</sup> *<sup>T</sup>*∞) <sup>=</sup> 0.094.

**Figure 9.** Comparison of temperature profile between WJ and the present work as a function of (**a**) time, and (**b**) tissue layer depth.

The comparison between WJ and the present work exhibited in Figure 9, for the temperature as a function of (a) time and (b) deep of the tissue layer, both fit in a good agreement. The comparison between the temperature distribution of the WJ model (40) and the present work (26) is shown in Figure 10. The heat source is applied on the central part of the tissue surface, as shown in Figure 2. The temperature distribution corresponds to a heat source application dimensionless time of *τ* = 0.09. The temperature distribution of the WJ model shows a higher evolution, indicating a higher heat transfer in the tissue. The difference between models according to temperature contours is less than 2%.

**Figure 10.** Comparison of temperature distribution between (**a**) WJ model, and (**b**) present work at dimensionless time. Just before the heat source was removed.

#### *4.3. Dynamical Test Simulations*

#### 4.3.1. One-Layer Tissue Analysis

To evaluate the different parameters involved in the ETC, a dynamic test is performed, which consists of applying a heat source on the tissue surface for a period of time, removing the source when the thermal relaxation process begins, the scheme and conditions are shown in Figure 2. This test does not reach the tissue degradation temperature, *T* < 44 ◦C [8,64].

Figure 11 presents the temperature profile as a function of time and tissue layer depth for one-layer (at *τ* = 0.09), by modifying the fractal coefficient *D<sup>T</sup>* in the range of 1–1.1, with porosity *φ* = 0.1 and Newtonian fluid *n* = 1. There are no significant changes when the heat source is applied, at the end of the relaxation process is when there is a difference between the WJ model and this work. This difference is due to the fact that there is a better energy transfer process in the present model, but there is no appreciable effect of *D<sup>T</sup>* due to the low tissue porosity.

**Figure 11.** Temperature profiles as a function of (**a**) time and (**b**) tissue layer depth, for different *D<sup>T</sup>* values at *τ* = 0.09.

Figure 12 shows the temperature profile as a function of time and tissue layer depth (at *τ* = 0.09) for one-layer tissue by modifying the porosity *φ* in the range of 0.01 to 0.5, with *D<sup>T</sup>* = 1.05 and Newtonian fluid *n* = 1. As the porosity increases, higher temperatures are reached when the heat source is applied, and lower temperatures are reached at the end of the relaxation process. This is due to increased perfusion in the tissue, on the other hand, blood has a higher thermal conductivity which improves heat transfer in the porous medium.

Figures 13 and 14 exhibit the temperature profile as a function of time and tissue layer depth (at *τ* = 0.09) one layer for two porosity values *φ* = 0.1 and *φ* = 0.5, respectively. Modifying the power-law index *n* in the range of 0.6 to 1 and *D<sup>T</sup>* = 1.05. In the case of porosity *φ* = 0.1, there are no significant changes when modifying the power-law index, due to low perfusion through the tissue. For porosity *φ* = 0.5, blood perfusion is increased and the effect of the power-law index is reflected both in the application of the heat source and in the thermal relaxation process. For a shear-thinning fluid *n* = 0.6, there is a greater dissipation of energy and therefore it reaches lower temperatures compared with a Newtonian fluid, in addition to experiencing a higher thermal relaxation due to a greater influence of the convective condition of the surface, as shown in Figure 14a. The effect of higher porosity is also reflected in the temperature profile as a function of tissue depth, as it presents distortions, mainly for the non-Newtonian fluid with *n* = 0.6, see Figure 14b.

**Figure 12.** Temperature profiles as a function of (**a**) time and (**b**) tissue layer depth, for different porosity *φ* values at *τ* = 0.09.

**Figure 13.** Temperature profiles as a function of (**a**) time and (**b**) tissue layer depth for different power-law index *n* at *τ* = 0.09 and porosity *φ* = 0.1.

**Figure 14.** Temperature profiles as a function of (**a**) time and (**b**) tissue layer depth, for different power-law index *n* at *τ* = 0.09 and porosity *φ* = 0.5.

4.3.2. Three-Layer Tissue Analysis

The structure of the skin is very complex and is generally considered in three layers epidermis, dermis, and hypodermis—as is shown in Figure 1. The thickness of these layers

varies depending on the location of the skin and many factors, such as age, sex, race, endocrine, and nutritional status of the individual [7–10].

Figure 15 shows the temperature profile as a function of time and tissue three-layer depth (at *τ* = 0.09), each layer with its own conductivity, *k<sup>e</sup>* = 0.25 W/m· ◦C, *k<sup>d</sup>* = 0.45 W/m· ◦C and *k<sup>h</sup>* = 0.2 W/m· ◦C, for epidermis, dermis, and hypodermis, respectively, [8]. The porosity *φ* varies uniformly throughout the tissue in the range of 0.01–0.5, with *D<sup>T</sup>* = 1.05 and non-Newtonian fluid *n* = 0.6. The shear-thinning effect increases the heat transfer in tissue, which is reflected in the temperatures reached when applying the heat source and in the thermal relaxation process. The effect of having three different conductivities can be seen clearly in the temperature profile as a function of tissue depth, where changes in the slope of the curve at the layer interfaces are presented, as is shown in Figure 15b. This behavior has been reported by Xu et al. [2].

**Figure 15.** Temperature profiles as a function of (**a**) time and (**b**) tissue three-layer depth with different porosities, *φ* = 0, 0.5, and 0.05 for epidermis, dermis, and hypodermis, respectively.

Figure 16 shows the temperature profile comparison between one- and three-layer tissue at different depths as a function of time and tissue depth, the properties of the three-layer tissue are those used in Figure 15. In the case of a one-layer tissue, the uniform porosity is *φ* = 0.5, for the case of three-layer tissue the following porosities are assigned *φ* = 0, 0.5, and 0.05 for epidermis, dermis, and hypodermis, respectively. These values were assigned according to the fact that the epidermis has very low porosity, the dermis contains the highest density of capillaries, and the hypodermis is composed of loose fatty connective tissue. The other variable values are *D<sup>T</sup>* = 1.05 and *n* = 0.6 [8]. Considering blood as a non-Newtonian shear-thinning fluid, in addition to different tissue properties, generates important modifications in the temperature profiles both in time and depth. In contrast, uniform properties throughout the biological tissue, as shown in Figure 16a,b.

Figure 16a shows red dots on temperature profiles used to indicate the time corresponding to the instantaneous temperature distribution presented in Figure 17. This figure shows the comparison of temperature distribution between the one- and three-layer tissues at different times. In the case of three-layer tissues, the effect of the different tissue properties on the temperature contours, which show small variations due to porosity depending on the area. In the case of one-layer tissues, the heat transfer is more similar to homogeneous solids, which have higher thermal conductivity at the top and bottom of the tissue compared with the three-layer model. Finally, the right half is overlapped on the left half of the temperature distribution for one- and three-layer tissue to demonstrate the asymmetry generated by the global effect of all the variables involved in the model, as shown in Figure 18.

**Figure 16.** Temperature profile comparison between one- and three-layer tissue as a function of (**a**) time at three different depths of the tissue and (**b**) depth of the three-layer tissue with different porosities *φ* = 0, 0.5, and 0.05 for epidermis, dermis, and hypodermis, respectively.

**Figure 17.** Comparison of the temperature distribution for one-layer tissue with *φ* = 0.5 and threelayer tissue with *φ* = 0, 0.5, and 0.05 for epidermis, dermis, and hypodermis, respectively. For different times from top to bottom (**a**) *τ* = 0.01, (**b**) *τ* = 0.09, (**c**) *τ* = 0.2, and (**d**) *τ* = 0.45. For a non-Newtonian fluid, *n* = 0.6.

**Figure 18.** The right half is overlapped on the left half of the temperature distribution for the oneand three-layer tissue. Temperature distribution was taken from Figure 17 at time *τ* = 0.09.

#### **5. Conclusions**

In this work, the ETC for human skin using the fractal scaling and Monte Carlo methods was obtained. These methods were used to describe the tissue as a porous medium, the blood was considered as Newtonian and non-Newtonian shear-thinning fluid. The numerical code developed was validated by comparing the ETC obtained by Weinbaum and Jiji [11] and the present work using one-layer tissue with the same porosity, straight capillaries, and Newtonian fluid. The difference between models according to temperature contours is less than 2%.

The ETC involves various parameters, such as fractal dimensions *D<sup>T</sup>* and *D<sup>f</sup>* , porosity, and the power-law index *n*; in order to evaluate these parameters, dynamical tests were performed, which consisted of applying a heat source on the tissue surface for a period of time—when the heat source was removed, the thermal relaxation process began. The effect of the main parameters on the temperature profiles as a function of time and tissue depth, for one- and three-layer tissue, besides temperature distribution, were presented. The main findings of this work are the following:


**Author Contributions:** Conceptualization, R.O.V. and G.R.-A.; methodology, R.O.V., G.R.-A., J.P.E. and R.M.-M.; software, R.O.V., G.R.-A., J.P.E. and A.R.-M.; validation, R.O.V. and G.R.-A.; formal analysis, R.O.V. and G.R.-A.; investigation, R.O.V., G.R.-A., J.P.E., A.R.-M. and R.M.-M.; resources, R.O.V., G.R.-A., J.P.E., A.R.-M. and R.M.-M.; data curation, G.R.-A.; writing—original draft preparation, R.O.V., G.R.-A. and J.P.E.; writing—review and editing, R.O.V., G.R.-A., J.P.E., A.R.-M. and R.M.-M.; visualization, R.O.V., G.R.-A. and J.P.E.; supervision, R.O.V.; project administration, R.O.V.; funding acquisition, R.O.V. and J.P.E. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors R.O.V. and J.P.E. would like to acknowledge the financial support SIP-20221783 and SIP-20221572.

**Acknowledgments:** Guillermo Rojas-Altamirano gratefully acknowledges support by the Instituto Politécnico Nacional of Mexico for the PhD scholarship.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**

