Abstract
A dynamic model predicting human thermal responses in cold, cool, neutral, warm, and hot environments is presented in a twopart study. This, the first paper, is concerned with aspects of the passive system:1) modeling the human body, 2) modeling heattransport mechanisms within the body and at its periphery, and 3) the numerical procedure. A paper in preparation will describe the active system and compare the model predictions with experimental data and the predictions by other models. Here, emphasis is given to a detailed modeling of the heat exchange with the environment: local variations of surface convection, directional radiation exchange, evaporation and moisture collection at the skin, and the nonuniformity of clothing ensembles. Other thermal effects are also modeled: the impact of activity level on work efficacy and the change of the effective radiant body area with posture. A stable and accurate hybrid numerical scheme was used to solve the set of differential equations. Predictions of the passive system model are compared with available analytic solutions for cylinders and spheres and show good agreement and stable numerical behavior even for large time steps.
 dynamic simulation
 human heat transfer
 asymmetric thermal environments
 exercise
 numerical modeling
in the past three decades, several multisegmental mathematical models of human thermoregulation have been developed (e.g. Refs. 17, 34, 41) and have become valuable tools in the hands of their authors in contributing to a deeper understanding of regulatory processes. These models, however, have not yet gained widespread use. Reasons for this may include lack of confidence in their predictive abilities, limited range of applicability, and poor modeling of the heat exchange with the environment. This paper and the paper in preparation make a contribution to research efforts to formulate a more precise, flexible, and universal model of the human thermoregulatory system.
The main stimulus for the model is the need to predict the thermal response of human beings to the conditions that they experience in and around buildings. Other wellestablished models can predict ambient air and surface temperatures, air speeds, and humidity; therefore, these parameters are the boundary conditions that will act on the thermoregulatory model.
The research sought to draw together the most appropriate models describing individual parts of the whole thermoregulatory system. These were then combined within a mathematical framework to produce a robust, computationally fast, and accurate model that will be of practical value in the design of buildings and urban environments. The research thus seeks to deliver the products of fundamental physiological research to scientists working in other disciplines.
From the mathematical point of view, the human organism can be separated into two interacting systems of thermoregulation: the controlling active system and the controlled passive system. The active system is simulated by means of cybernetic models predicting regulatory responses, i.e., shivering, vasomotion, and sweating, as discussed elsewhere (16). The passive system, the subject of this paper, is modeled by simulating the physical human body and the heattransfer phenomena occurring in it and at its surface. Comparisons between the model's predictions and known analytic solutions are made. Only the complete model, i.e., the passive system plus active system, can be compared with field measurements of human responses to environmental conditions and with the predictions of whole body models produced by others; such comparisons are, therefore, a subject of the upcoming paper.
The metabolic heat is produced within the body. This heat is distributed over body regions by blood circulation and is carried by conduction to the body surface, where, insulated by clothing, the heat is lost to the surroundings by convection, radiation, and evaporation (see Fig. 1). Additionally, some heat is lost to the surroundings by respiration. To simulate adequately all these heattransport phenomena, the present model of the passive system accounts for the geometric and anatomic characteristics of the human body and considers the thermophysical and the basal physiological properties of tissue materials.
Whereas the model was inspired by the desire to evaluate environmental conditions in buildings, it also has applications in other areas, e.g., in the car industry, in textile industries, in the aerospace industry, in meteorology, in medicine, and in military applications. In these disciplines, the model can serve for research into human performance, thermal acceptability and temperature sensation, clinical and therapeutic treatments, safety limits, etc.
BODY CONSTRUCTION
The passive system of the model was developed by using thermophysical and basal physiological properties of tissue materials, based on literature (10, 17, 40). The humanoid was formulated to represent an average human (see Table 1 for overall body data). The body was idealized as 15 spherical or cylindrical body elements: head, face, neck, shoulders, arms, hands, thorax, abdomen, legs, and feet. The use of lumped data (e.g., if body elements are divided into a core and a shell) is prone to errors during temperature transients and was thus avoided. In accordance with former studies (17), a division was therefore made whenever a significant change of body tissue properties occurred. As a result, the present multilayer model consists of annular concentric tissue layers and uses seven different tissue materials: brain, lung, bone, muscle, viscera, fat, and skin. Each of these is subdivided into one or more tissue nodes (see Table 2).
In accordance with the work of Weinbaum et al. (39), the skin was modeled as two layers with distinctly different physiological properties: the inner skin and the outer skin. The inner skin is ∼1 mm thick and simulates the cutaneous plexus, a region in which metabolic heat is generated and blood is perfused (see Table 2). The outer layer of skin is of similar thickness; however, it contains neither any heat source nor any thermally significant blood vessels. This superficial cutaneous layer plays a role when evaporative heat losses through the skin are modeled. The layer contains sweat glands and simulates the vapor barrier for moisture diffusion through the skin.
Asymmetric removal of bodily heat is often experienced in practice. It may be due to inhomogenous ambient conditions, e.g., solar irradiation, cold surfaces, draught effects, nonuniform clothing, or to physiological effects on the skin, such as local differences in evaporative heat losses, etc. Except for the face and shoulders, each body element was, therefore, subdivided spatially into sectors. Most of the body elements were divided into three sectors: anterior, posterior, and inferior (see Table 2). Anterior and posterior segments permit the treatment of lateral environmental asymmetries. Inferior segments account for body sides, which are “hidden” by other body parts (e.g., section AA′ in Fig. 1) and thus have reduced radiant heat exchange with the environment.
The sectors were coupled thermally by introducing a core element around the cylinder axis (see Fig. 1, section AA′) or the midpoint of the sphere. In the model, except for the head, the radius of the core was defined to be identical to the radius of the innermost tissue material layer in each body element. In extremities, for instance, the radius of this domain is identical to the radius of the bone. However, the brain, which has an outer radius of r = 8.6 cm (Table 2), was given a core of 4cm radius.
Glossary
Symbols
γ, δ, ζ
 α
 Shortwave absorptivity
 β
 Factor of the blood perfusion rate (W ⋅ K^{−1} ⋅ m^{−3})
 Δ
 Difference
 ς
 StefanBoltzmann constant (W ⋅ m^{−2} ⋅ K^{−4}) (=5.67051 × 10^{−8}W ⋅ m^{−2} ⋅ K^{−4})
 ε
 Longwave emissivity
 γ, δ, ζ
 Timeindependent coefficients arising from the numerical formulation of the bioheat equation
 λ_{H2}_{O}
 Heat of vaporization of water (= 2,256 × 10^{3}J/kg^{−3})
 μ_{bl}
 Proportionality factor (K^{−1})
 η
 Meachanical efficiency of the human body
 ρ
 Density (kg/m^{3})
 ψ
 View factor
 Ω
 Linear operator
 ϕ
 Angle of a body element sector (°)
 Σ
 Sum
 κ
 Geometry factor of an isothermal core domain of a body element
 ω
 Geometry factor of the bioheat equation
 ∂
 Partial derivative
 A
 Area (m^{2})
 a
 Distribution or weighting coefficient
 act
 Activity level (met)
 b
 Regression coefficient
 B_{x}
 Arterial blood temperature parameter arising from h_{x}
 C
 Heat loss by convection (W)
 c
 Heat capacitance (J ⋅ kg^{−1} ⋅ K^{−1})
 Cc
 Coefficient of the blood pool temperature column of the major coefficient matrix
 CCX
 Countercurrent heat exchange
 CO
 Cardiac output (l/m)
 CplC
 Coupling coefficient of the major coefficient matrix
 E
 Heat loss by evaporation (W)
 f_{cl}
 Overall clothing surface area factor
 f*_{cl}
 Local value of f_{cl}
 f_{eff}
 Effective radiant area of the body
 h
 Surface heat transfer coefficient (W ⋅ m^{−2} ⋅ K^{−1})
 h_{x}
 Countercurrent heat exchange coefficient (W/K)
 H
 Internal whole body workload (W)
 i_{cl}
 Overall moisture permeability index from the skin to the clothing surface
 i*_{cl}
 Local value of i_{cl}
 I_{cl}
 Overall intrinsic insulation from the skin to the clothing surface (m^{2} ⋅ K ⋅ W^{−1})
 I*_{cl}
 Local value of I_{cl}(m^{2} ⋅ K ⋅ W^{−1})
 k
 Conductivity (W ⋅ m^{−1} ⋅ K^{−1})
 L
 Lewis constant (L_{a} = 0.0165) (K/Pa)
 Lc
 Coefficient of the blood pool temperature line of the major coefficient matrix
 m
 Mass (kg)
 M
 Whole body metabolism (W)
 MRT
 Mean radiant temperature (°C)
 N
 Number of nodes
 P
 Water vapor pressure (Pa)
 q
 Heat flux density (W/m^{2})
 q_{m}
 Metabolic rate (W/m^{3})
 Q
 Heat flux (W)
 r
 Radius (m)
 R_{e}
 Moisture permeability resistance (m^{2} ⋅ Pa ⋅ W^{−1})
 rh
 Relative humidity (%)
 s
 Shortwave irradiation (W/m^{2})
 t
 Time (s)
 T
 Temperature (°C)
 T*
 absolute temperature (K)
 U*
 Local clothing insulation coefficient from the skin to the ambient air (W ⋅ m^{−2} ⋅ K^{−1}) (or) (W ⋅ m^{−2} ⋅ Pa^{−1})
 v
 Airspeed (m/s)
 V
 Volume (m^{3})
 w
 Perfusion rate (s^{−1})
 wt
 Wetted area ratio
Subscripts and superscripts
impl
 −
 Last node next to an interface between two inhomogenous tissue shells
 +
 First node next to an interface between two inhomogenous tissue shells
 0
 Thermal neutrality
 air
 Air
 acc
 Accumulation
 bas
 Basal value
 bd
 Body
 bl
 Blood
 bl,a
 Arterial blood
 bl,p
 Central blood pool
 bl,v
 Venous blood
 c
 Convection
 cl
 Clothing
 CN
 CrankNicolson
 Du
 Dubois
 e
 Evaporation
 eff
 Effective
 expl
 Explicit
 f
 Fabric
 frc
 Forced
 ifc
 Interface
 j
 Running number variable
 hy
 Hypothalamus
 impl
 Implicit
 m
 Mean value
 mix
 Mixed
 mus
 Muscle
 n
 Total number
 nat
 Natural
 o
 Operative
 osk
 Outer skin layer
 r
 Node number
 R
 Long wave radiation
 re
 Rectal
 rsp
 Respiration
 sat
 Saturated
 sf
 Surface of body sector
 sh
 Shivering
 sk
 Skin
 sR
 Shortwave radiation
 sr
 Surrounding
 sw
 Sweating
 t
 Time step number
 ty
 Tympanic
 w
 Work
 x
 Countercurrent heat exchange
HEAT TRANSFER WITHIN THE TISSUE
The heat transport mechanisms occurring in the living tissue have been formulated by Pennes (28) in the socalled “bioheat equation.” This differential equation describes the heat dissipation in a homogenous, infinite tissue volume
According to Eq. 1, from left to right, the radial heat flow from warmer to colder tissue regions [heatconduction term, wherek is the tissue conductivity (W ⋅ m^{−1} ⋅ K^{−1}), T is tissue temperature (°C), r is radius (m), ω is a geometry factor (dimensionless); ω = 1 for polar coordinates and ω = 2 for spherical coordinates (head)] is overlayed by metabolismq _{m}(W/m^{3}) and blood perfusion [heatconvection term, where ρ_{bl} is density of blood (kg/m^{3}), w _{bl} is blood perfusion rate (s^{−1}), c_{bl} is heat capacitance of blood (J ⋅ kg^{−1} ⋅ K^{−1}), and T_{bl,a} is arterial blood temperature (°C)]. This combined effect is balanced by the storage of heat within the tissue mass [righthand side of Eq. 1, where ρ is tissue density (kg/m^{3}), c is tissue heat capacitance (J ⋅ kg^{−1} ⋅ K^{−1}) and t is time (s)].
The bioheat equation was applied to all tissue nodes by using the appropriate material constants k, ρ, and c, the basal heatgeneration term q _{m}, and basal blood perfusion rate w _{bl} for each tissue layer (Table 2). It should be noted that the magnitudes of the physiological variablesq _{m} and w _{bl} are affected by responses of the active system, as described elsewhere (16). The arterial blood temperature T_{bl,a} in Eq. 1 arises from the actual overall thermal state of the body and is obtained by simulating the human blood circulatory system.
Heat conduction.
The heatconduction term in Eq. 1 considers temperature variations in the radial direction. Angular heat flows were also neglected in the present model, even though the body elements were divided into sectors, and environmental asymmetries were considered. There were two reasons for this simplification. First, at the periphery (i.e., predominantly in the skin), although significant angular temperature gradients may occur due to asymmetric boundary conditions, the total amount of the heat conducted angularly from one skin sector to another is negligible, compared with the total amount of heat transferred radially through a sector band. This is because peripheral sectors have large shell areas, compared with the small interface areas between adjacent skin sectors. Second, ambient influences are insignificant in internal tissue layers where heat dissipation is dominated by metabolism and by blood circulation (28,33). Therefore, inhomogeneities of q _{m} andw _{bl} between individual organs of the abdominal viscera (40), for instance, which can hardly be considered by symmetrical geometry models, govern local temperature profiles within the body and swamp the effect of environmental asymmetries.
Metabolism.
In the model, the metabolic heat production of a finite tissue volume is treated as a sum of the basal value q
_{m,bas,0} and the additional heat Δq
_{m}, which may be produced by local autonomic thermoregulation, while the subjects are exercising and/or shivering
The basal metabolic rates of individual tissue materialsq _{m,bas,0} were obtained from literature (10, 17, 40) and are listed in Table 2. In thermal neutrality (T_{air} = 30°C; for a detailed specification see Table3), where no regulation occurs, the overall basal metabolism M _{bas,0} = ∫q _{m,bas,0} dV of the passive system results in a value of 87 W, which agrees with the (standardized) whole body metabolism of a reclining average human, according to ASHRAE Standards 55 (1).
In muscles, the extra heat Δq
_{m} may contain three components; i.e., changes in the basal metabolism (Δq
_{m,bas}) and in additional metabolism by shivering and working (q
_{m,sh} andq
_{m,w}, respectively)
The change in the basal metabolism Δq
_{m,bas} is the difference between the actual basal rate and the basal rate corresponding to neutral thermal conditions. It appears also in nonmuscular tissues where q
_{m,bas,0} is present and where the tissue temperature differs from its setpoint T_{0}. This local regulation arises from the van't Hoff Q_{10}effect with a sensitivity coefficient of 2 (40), which reflects the dependence of biochemical reactions on local tissue temperature T
The workload H in Eq. 5
represents that part of the overall extra energy which is produced by muscular activity but which does not appear as external work. It is obtained as the difference between the actual overall heatgeneration rate remaining in the body and the basal metabolism M
_{bas,0}
The calculation of the actual metabolic heat production utilizes the activity level units act (measured in met) as the rate between the body heat produced at that activity and the average metabolism when the subject is seated, defined as 58.2 W/m^{2} (2), i.e., 1 met. In Eq. 6, the whole body metabolism is calculated from the basal metabolic rate M _{bas,0} for a reclining, resting individual, corresponding to act_{bas} = 0.8 met. The act units, used as (timedependent) input variables to the model, are commonly known and can easily be obtained by the users from standard literature (e.g., Refs. 1, 2).
The term in parentheses (1 − η) in Eq. 6
represents that fraction of the overall metabolic heat actually produced which remains for warming the body, i.e., which is not converted into external mechanical power. The human mechanical efficiency η is not constant in reality but rises with increasing activity levels. For the present model, we developed the following formula for η, by means of regression analysis from measurements of Wyndham et al. (42)
The regression coefficients were b _{1} = 0.39 ± 0.13 met^{−1} and b _{0} = −0.60 ± 0.28, with a correlation coefficient of r = 0.86. This equation reflects the behavior of the human work efficacy during normal activities [e.g., stepping or cycling (42)]. At low activity levels (act < 1.6 met), where humans work most inefficiently, η is zero. At higher activity levels (∼1.6 < act < 5.0 met), η rises almost linearly with increasing activity level. At high activity levels, η begins to approach a theoretical maximum of ∼0.2; i.e., ∼20% of the whole body metabolism is transferred into external work.
Blood circulation.
Blood circulation is of vital importance for the dissipation of heat within the human body. The dissipation effect was examined by simulations using the present model without blood circulation: the core temperature of the brain, which has a high metabolic rate (Table 2), exceeded 73°C when the passive system was exposed to steady neutral ambient temperatures of 30°C.
In the present model, the human blood circulatory system was simulated as three main components: 1) the central blood pool, 2) countercurrent heat exchanges CCX, and 3) pathways to individual tissue nodes. Body elements are supplied with warm blood from the central pool by the major arteries. However, before the tissue of the extremities and shoulders are perfused, the blood is cooled by heat lost to the countercurrent bloodstreams in the adjacent veins. The arterial blood exchanges heat by convection in the capillary beds, according to the bioheat Eq. 1. The venous blood then collects in the major veins and is rewarmed by heat from the adjacent arteries as it flows back to the central blood pool. It mixes with blood from other elements to produce the new central blood pool temperature.
The blood perfusion term of the bioheat Eq. 1 is based on Fick's first principle, assuming that heat is exchanged with the tissue in the capillary bed and that no heat storage occurs in the bloodstream. The blood perfusion term has been critically discussed in the literature (see, e.g., Ref. 11) and has been subjected to various improvements (13, 33). Also, alternative models of microvascular structures have been developed (see, e.g., Ref. 39). However, the traditional blood perfusion term seems preferable in whole body models. It is simple and produces an accuracy that is comparable (9) to more detailed models that require comprehensive data on the vascular architecture of tissues.
The simulation of the countercurrent heat exchange CCX was introduced into the model in an attempt to obtain a more realistic distribution of the arterial blood temperature T_{bl,a}, instead of assuming a constant arterial blood temperature for all body elements equal to the temperature of the central blood pool T_{bl,p}. A number of models of the CCX between pairs of adjacent vessels exist in the literature [for a review see, e.g., Eberhart (13)], but these are seldom incorporated into whole body models.
Assuming mass continuity in blood vessels, the decrease in the temperature of the blood in the arteries of an element T_{bl,p} − T_{bl,a} corresponds to the increase in temperature T_{bl,v,x} − T_{bl,v} of the blood in the veins of the same body element after CCX
The computation of the central blood pool temperature T_{bl,p}utilized a special calculation process (see numerical methods). Because the bioheat equation assumes that capillary blood reaches an equilibrium with the surrounding tissue, T_{bl,v} of a body element in Eq. 10
is obtained as follows
The heatexchange coefficients h_{x} in Eq. 10 are zero for central body elements (17) (i.e., no countercurrent heat exchange). Few reliable data on h_{x} of extremities exist in literature. Therefore, h_{x} values of extremities and shoulders were estimated by a trialanderror procedure, in which the coefficients were adjusted to obtain agreement between local skin temperatures as predicted and measured. For legs, feet, arms, and hands, h_{x} were obtained from measurements by Mayer and Schwab (21) and Nielsen and Pedersen (25). An experiment of Budd and Warhaft (6) provided data for the shoulders. The individual coefficients h_{x} are listed in Table 2.
In thermal neutrality, tissues are supplied at basal perfusion ratesw
_{bl,0} (see Table 2). In nonneutral conditions or during work, the blood flows, w
_{bl} in Eq. 1,vary with changes in regional metabolic rates q
_{m}. Experimental estimations (30) indicate a linear relationship betweenw
_{bl} and q
_{m}. Instead of dealing with blood perfusion rates w
_{bl}, the present model considers directly their energy equivalent β = ρ_{bl}c_{bl}
w
_{bl}(W ⋅ m^{−3} ⋅ K^{−1}). The increased demand for oxygen is therefore accounted for by calculating the change in the gain factor Δβ = ρ_{bl}c_{bl}Δw
_{bl} as a function of the change in metabolism Δq
_{m} using a proportionality constant μ_{bl} = 0.932 (K^{−1}), which was obtained from Stolwijk (34)
The very variable blood perfusion rates in the cutaneous plexus arise from human temperatureregulation mechanisms, as discussed elsewhere (16).
HEAT EXCHANGE WITH THE ENVIRONMENT (BOUNDARY CONDITIONS)
At the body surface, heat is exchanged by convection with the ambient air q
_{c}, by radiation with surrounding surfacesq
_{R}, by irradiation from hightemperature sourcesq
_{sR}, and by evaporation of moisture from the skinq
_{e}. The rate of heat exchange varies over the body and is affected by the clothing ensemble worn. Hence, in the model, heat balances were established for each skin sector of each body element as boundary conditions to Eq. 1
. In general, the net heat flux q
_{sk} (W/m^{2}) passing the surface of a peripheral sector is equivalent to the sum of these individual heat exchanges
The elements of this equation are considered in detail in following sections.
Convection.
The convective heat exchange q
_{c} between a skin sector of surface temperature T_{sf} and ambient air of temperature T_{air} was modeled by considering both natural and forced convection using combined convection coefficients h_{c,mix}
In the present model, the convection coefficients h_{c,mix}(W ⋅ m^{−2} ⋅ K^{−1}) are a function of the location on the body, the temperature difference between the surface and the air, and the effective airspeedv
_{air,eff} (m/s)
The basis for h_{c,mix} are the experiments of Wang (3638), who measured local convective heat losses on a heated fullscale manikin with a realistic skin temperature distribution for both natural and forced convection. Coefficientsa
_{nat}, a
_{frc}, anda
_{mix} (see Table 2) were obtained by means of regression analysis in which Wang's results were modified and transformed into expressions more appropriate for computer simulation, as shown by Eq. 17. To check the reliability of the convective heatexchange model, the overall, i.e., the mean, convection coefficient of the human body, h_{c,m}(W ⋅ m^{−2} ⋅ K^{−1}), was obtained by integrating the local convective heat lossesq
_{c} = h_{c,mix}(T_{sk} − T_{air}) over the (naked) body surface A
_{Du}
The overall coefficient was computed both for conditions of the natural convection (v _{air,eff} < 0.05 m/s) and for forced convection. Figure 2 compares the simulation results for the free convection with the data of Fanger (14) and Bach (3) and the data quoted by McIntyre (24), i.e., free convection coefficients for a 2mhigh cylinder obtained from Birkebak (5) and coefficients obtained by Rapp (29), which are recommended for sedentary people by McIntyre (24). Simulated forced convection coefficients were compared with formulas used by Fanger (14), Bach (3), Colin and Houdas (12), and Nishi (26) and are illustrated in Fig.3. As apparent from Fig. 3, the model's overall coefficient for forced convection correlates better with results of more recent literature (3) indicating higher values (see also Refs. 8 and 21), compared with earlier investigations.
The effective air velocity v _{air,eff} in Eq.17 is a composite of the local room airspeed and the relative airspeed due to body motion. Earlier investigations of walking subjects (26) suggested an increase in h_{c,m} with increasing body movement. More recent experiments (8), however, have illustrated constant or even decreasing h_{c} coefficients during walking. Similarly, simulations using the present model did not show evidence of functional dependency of the convective heat loss on the activity level. Obviously, h_{c} is affected by the type of activity rather then by its level. Hence, instead of introducing a fixed relationship between h_{c} and the activity level, which could produce systematic errors in predictions, users of the model may modify the input variable v _{air,eff} appropriately, if detailed information on h_{c} for a specific task is available.
Radiation.
The exchange of heat by longwave radiation q _{R} is usually of similar importance in the heat balance of the human body as the heat exchange by convection.
In asymmetric environments, q
_{R} of a body element sector represents the sum of the partial heat exchanges between this sector and the surrounding structures (walls, windows, etc.), which may have different surface temperatures. Thus the exact simulation ofq
_{R} requires the computation of view factors (which describe the influence of geometry on the radiative heat exchange) between the surface sectors of the body and each surface segment of the surroundings. Unfortunately, calculation of the view factors between the curved surfaces of the body and surrounding structures is very computer intensive. Hence, the concept of directiondependent average temperatures for the surrounding structures was adopted. In the model, either T_{sr,m}, the mean temperature of the surrounding surfaces, or mean radiant temperature MRT can be used. Each of these will differ from one body sector to another, and they may vary with time. The value of T_{sr,m} is defined as the temperature of a fictitious uniform envelope “seen” by a sector, which causes the same radiative heat exchange with the sector as the actual asymmetric enclosure. MRT is defined similarly, but it refers to a fictitious uniform black envelope. Employing T_{sr,m}, and rearranging the StefanBolzmann law by introducing the (local) radiative heatexchange coefficient h_{R}(W ⋅ m^{−2} ⋅ K^{−1}), the energy exchange is calculated by
The longwave emissivity of the body segments ε_{sf} depends on the covering material. The values for skin and hair in Table 2 were obtained from Refs. 2 and 17. Clothing usually has a value close to 0.95 (14). The emissivity of the surroundings, which is an input variable into the model, is typically ε_{sr} = 0.93 for indoor spaces. Note that if MRT is used in place of T_{sr,m}, ε_{sr} has to be set to unity (black envelope).
The concept of using mean surrounding temperatures T_{sr,m}makes it possible to introduce the view factors ψ_{sfsr}between sectors and the surroundings seen by them. The view factors ψ_{sfsr} of the cylindrical and spherical body elements (Table 2) were calculated—for each sector separately—by a specially developed finiteelement program (15). They vary from 0.1 to unity, depending on the degree to which individual skin sectors are “hidden” by other body parts. View factors are given for postures with stretched limbs (reclining, standing activities) and for the seated position that is characterized by a decreased bodyradiation area. As a check on the reliability of the view factors, the effective radiant area ratio of the (nude) body f _{eff} = ∫ ψ_{sfsr}dA _{sk}/A _{Du}was computed. The resultant whole body values off _{eff} = 0.80 for the standing position andf _{eff} = 0.74 for the seated position are somewhat higher than provided, e.g., by Fanger (14), but they are in excellent agreement with results of recent, more detailed investigations (19).
The mean surrounding temperatures are provided together with other environmental variables in a socalled climateinput file. For inhomogeneous environments, T_{sr,m} should be given separately for the half of the room in front of the person (to be applied to the anterior sectors of the body elements) and the half behind the person (for posterior sectors). However, sometimes, it might be useful to determine T_{sr,m} for additional directions. For example, T_{sr,m} would be calculated separately for the upper part of the room when investigating the radiative effects of heating/cooling ceilings, which may cause occupants to perceive local discomfort to the head and shoulders.
While the MRT of a half room can be measured directly, the corresponding T_{sr,m} must be calculated from measurements of individual surface temperatures. A simple approach to calculate T_{sr,m} is to weight the surface temperatures T_{sr, j} by the corresponding surface areasA
_{sr, j}
A more detailed method is to consider the effect of the body's location in the room by using the view factors ψ_{bdsr, j}
The overall view factors ψ_{bdsr,j} between the body and wall segments have been measured experimentally by means of photographic projection techniques for both standing and sitting postures (14) and can either be taken from graphs or calculated by using appropriate algorithms (7) and computer models (35).
The irradiation of the body by hightemperature sources (sun, fireplaces, etc.) was also considered in the formulation of the passive system. The shortwave radiation intensity and its direction is a characteristic of the environment and is provided, together with other environmental variables, in the climateinput file. The amount of heatq
_{sR} (W/m^{2}) absorbed at a sector surface is taken into account in the heat balance of superficial body element sectors by the term
Clothing insulation.
Because clothing plays an important role in the behavioral thermoregulation of human beings, efforts were made to model garments in some detail. Unfortunately, standard literature provides only overall clothing parameters, i.e., the intrinsic clothing insulation I_{cl}, the clothing area factor f _{cl}, and the moisture permeability index i _{cl} (for exact definitions see, e.g., Ref. 22). These do not reflect the real insulating properties of the garments but provide characteristics of a uniform (imaginary) clothing layer covering the entire body. Therefore, a simple computer model was developed that uses the present passive system and resimulates the experimental procedures for measuring the overall insulation characteristics of garments undertaken by McCullough et al. (22, 23). The program computes the required local insulation values I*_{cl},f *_{cl}, andi*_{cl} of trousers, shirts, underwear, socks, etc., from measured uniform layer data I_{cl},f _{cl}, and i _{cl} for these garments and from the evaporative resistance of the fabric R_{e,f}. These local values include the insulation effect provided by any layer of air trapped between the skin and the fabric. A database on clothing items was created, which was used to compose clothing ensembles in the present multisegmental model.
The local effective heattransfer coefficientU*_{cl}(W ⋅ m^{−2} ⋅ K^{−1}) of a jlayered clothing ensemble worn on a bodyelement sector is computed in the model as
The corresponding evaporative coefficientU*_{e,cl}(W ⋅ m^{−2} ⋅ Pa^{−1}) is obtained by using the local values i*_{cl}and I*_{cl} of individual clothing layers applied to the body sector, f *_{cl} of the outer clothing layer, the local convection coefficient h_{c,mix}, and the Lewis constant for air L
_{air} = 0.0165 K/Pa (24)
As a check on the validity of the new clothing model, the overall insulation I_{cl} of clothing ensembles and of individual garments was computed by using a similar procedure to that employed for obtaining the overall convection coefficient h_{c,m}. However, ∫q _{c}dA _{sk},A _{Du}, and T_{air} in Eq. 18 were replaced by the integral of the local heat fluxes through the clothing ∫ q _{cl}dA _{cl}, the overall surface area A _{cl} of the clothed body, and the corresponding mean surface temperature T_{cl,m}, respectively. The comparisons with results of experiments on clothed subjects showed good agreement with published data. For instance, simulating the thermal comfort experiment of Olesen and Fanger (27), in which the subjects wore the socalled “KSUuniform,” I_{cl}= 0.6 clo, i _{cl} = 0.34, the overall insulation was predicted to be I_{cl} = 0.59 clo andi _{cl} = 0.34.
Evaporation.
The present skin evaporation model ensures heat and mass transfer balances at each skin sector of each body element. The latent energy transport from a skin sector of area A
_{sk} is given by
The lefthand side of Eq. 26 represents the net energy transfer as driven by the evaporative potential between skin and air, where P_{sk} is the water vapor pressure at the skin surface, P_{air} is the vapor pressure of the ambient air, andU*_{e,cl} is the resultant evaporation coefficient of garments covering the sector.
The first righthand term of Eq. 26 considers the evaporation of sweat from the skin surface, whereby λ_{H2} _{O} = 2,256 kJ/kg is the heat of vaporization of water and dm_{sw}/dt is the rate of sweat production over A _{sk}, as elicited by the active system (16).
The last term of Eq. 26
describes the heat transport by moisture diffusion through the skin where P_{osk,sat} is the saturated vapor pressure within the outer skin layer of the body sector considered. In this superficial cutaneous layer, where sweat glands are located, moisture is always present. Hence, the vapor pressure at this location is set to its saturated value P_{osk,sat} depending on the tissue temperature of the outer skin, T_{osk} (24)
The net evaporative heat loss of a skin sector arises from the actual vapor pressure P_{sk} (Pa) found at the skin surface. This quantity is calculated by using a rearranged Eq. 26
The present evaporation model accounts for the storage of sweat liquid at the skin surface, as described by Jones and Ogawa (20). Moisture is accumulated when P_{sk}, as computed from Eq. 28,exceeds the saturated vapor pressure P_{sk,sat}, which is calculated from T_{sk} by using an equation similar to Eq.27. The rate of moisture storage dm_{acc}/dt(kg/s) is then equivalent to the rate of moisture production dm_{sw}/dt minus the rate of moisture evaporation
According to the experimental results quoted by Jones and Ogawa (20), the maximum amount of sweat liquid on the skin is limited to about m_{acc}/A _{sk} = 35 g/m^{2}. Quantities exceeding this threshold run off and are not considered in the model.
Respiratory heat losses.
Most heat is lost through the body surface; however, heat exchange with the environment also occurs by respiration. The model for respiratory heat loss was obtained from the work of Fanger (14), taking into account the loss by evaporation and convection. The latent heat exchange is calculated according to Ref. 14 from the pulmonary ventilation as a function of the whole body metabolism [i.e., ∫q
_{m}dV (W)]; latent heat of vaporization of water; and the difference between humidity ratio of the expired and inspired air, which depends on the ambient air temperature T_{air}(°C) and the partial vapor pressure of the ambient air, P_{air} (Pa)
Because the enthalpy of the expired air still depends to a certain degree on the condition of the inspired air, E _{rsp}appears as a function of T_{air} and P_{air}, allowing E _{rsp} to be calculated for a wide range of ambient conditions.
The dry heat loss of respiration due to the temperature difference between expired and inspired air can be expressed (14) as a function of the pulmonary ventilation rate and the temperature and vapor pressure of the ambient air, T_{air} and P_{air}, respectively
The total respiratory heat loss E
_{rsp} + C_{rsp} was distributed over body elements of the pulmonary tract: the lung and the muscle layers of the face and of the neck. The distribution coefficients a
_{rsp} were estimated from the work of Scherer and Hanna (31). This approach recognizes that inspired air is mainly conditioned in the nasal cavity where it reaches ∼60% of its final enthalpy deep in the lung. By considering also the rewarming of the nasal region by expiration, the following distribution percentages were obtained: 25% in the outer face muscles, 20% in the inner face muscles, 25% in the muscle band of the neck, and 30% in the lung. These heat losses appear in the nodal heat balances for these tissues as the volume derivative in the heatgeneration termq
_{m} in Eq. 1
NUMERICAL METHODS
In the model, each sector of each tissue shell was divided into nodes. These were unequally spaced in the radial direction, being closer together in the outer regions where temperature gradients are generally the steepest and where, therefore, the largest numerical errors could occur (17). The number of nodes per tissue shell N can be chosen arbitrarily in the model; it is only limited by the capacity of the computer memory. The current nodal configuration is shown in Table2.
A finitedifference scheme was used to discretize Eq 1.The partial derivatives with radius were approximated by the central difference method, which provides improved (i.e., secondorder) accuracy, compared with the first order of the forward or backward difference methods. The partial derivative for the temporal change was approximated by using a hybrid scheme, the CrankNicolson method, which also has a second order of accuracy. This method remains unconditionally stable, so the model can tolerate time steps of variable length. This means that the model can be more easily integrated with other dynamic simulation programs (e.g., of the building or vehicle environment).
The resulting discrete formulation of the bioheat equation was equally applicable to the cylindrical and to the spherical body parts. The passive system was organized as a structured system of timeindependent “conduction” matrices, which collected all thermophysical tissue constants, and timedependent “blood” matrices, the coefficients of which had to be computed for each time step of the program. Details are given in the . A quick solver was developed, which worked efficiently because it considered only the nonzero cells in the matrices.
RESULTS
Because the passive system does not contain models for the thermoregulatory feedback mechanisms that exist in human beings, it was not possible to validate the model by comparisons with field measurements. It was, however, possible to verify the numerical accuracy of the passivesystem model by comparing its predictions with existing analytic solutions for heat flow in homogeneous cylinders and spheres. Of most concern was the accuracy of the predicted temperatures, but the correct formulation of the boundary conditions was also examined. Tests were performed for a variety of thermal conditions in which the passive system was exposed to cold air. Because of producing steep temperature gradients in the tissues, these conditions reveal computational errors better than the warmair or hotair exposures (17).
Analytic solutions for conduction in cylinders and spheres as well as for conduction superimposed by internal heat generation (18) were employed to verify the model in steadystate conditions. These tests indicated that the model predicts steadystate tissue temperatures, which are so close to the known analytic results that the mean error is typically <0.02 K.
The performance of the model in transient conditions was of special interest. An analytic solution for transient heat conduction (without heat generation and fluid perfusion) in spheres derived by Gröber et al. (18) was compared with the predictions of the model for the head sphere made of brain tissue and setting q _{m} = β = 0. At t = 0, the sphere, which had an initially homogenous temperature of 37°C, was exposed to ambient temperatures of T_{air} = 0°C. Because there was no internal heat source, the sphere continuously cooled down. Good agreement with the solution of Gröber et al. was obtained for all temperatures in the sphere, as shown in Fig. 4 for three nodes (core, periphery, and a node in between). An error analysis of all nodal temperatures and times indicated a mean error of ‖ΔT‖ = 0.19 K, with the largest deviations from analytic results arising in the node next to the core, ‖ΔT‖ = 0.45 K (core: ‖ΔT‖ = 0.15 K), and the smallest errors occurring at the periphery, ‖ΔT‖ = 0.04 K. Simulations with a time step of 60 s did not significantly increase the accuracy for the given nodal configuration (see Table 2). The numerical scheme employed provided reasonable accuracies for both polar and spherical coordinates (mean error of ‖ΔT‖ < 0.5 K), even for time steps (Δt) of 1 h.
In another severe test of the model, the calculated temperature within a homogenous cylinder of muscle tissue was compared with an analytic solution of the bioheat equation quoted by Eberhart (13). The cylinder is initially in thermal equilibrium with the environment at T = T_{air} = 0°C. At t = 0, the cylinder is suddenly supplied with a constant rate of blood flow (w _{bl} = 0.001 s^{−1}, T_{bl,a} = 37°C) and source of metabolic heat (q _{m} = 600 W/m^{3}). The temperature of the cylinder therefore rises quickly before converging asymptotically to the final steady value. The predicted temperatures for three nodes in the cylinder, the core, the periphery, and a node between them, are compared with the analytically derived values in Fig.5. The mean error for all temperatures at Δt = 5 min was ‖ΔT‖ = 0.04 K. The reverse error behavior to the previous example was observed: larger deviations from the exact values occurred at the surface (‖ΔT‖ = 0.11 K) and the lower mean errors occurred in the core (‖ΔT‖ = 0.04 K).
The effect of timestep length was investigated by using the previous test but with a higher blood perfusion rate (w _{bl} = 0.005 s^{−1}), to ensure an even more rapid heating of the cylinder (reaching ∼95% of the steadystate value after 10 min). The simulation was undertaken with a time step of 10 min, and the predictions were compared with the exact solution (Fig. 6). There was an initial “overshoot” of the predicted temperature in the outermost cylinder node, but the temperature oscillations were strongly damped, and the predicted values quickly approached the analytic solution. This demonstrates the stability of the model even for disproportionally large time steps.
Finally, the temperature distribution within, and at the surface of, the whole passive system was calculated for conditions of “thermal neutrality” of T_{air} = 30°C (4) (for a detailed specification see Table 3). The simulation of the nonregulated, reclining, and unclothed human body provided a mean skin temperature of T_{sk,m} = 34.4°C, and a head core temperature (hypothalamus) of T_{hy} = 37.0°C, values that agree with published data of T_{sk,m} = 34.4°C (14) and T_{ty} = 37.0°C (4) for reclining subjects exposed to thermoneutral conditions. Some further simulation results of this steadystate exposure are listed in Table4.
Conclusions.
The heat exchange with the environment plays a key role in the thermal state of the human body. Development of models for the passive system of human thermoregulation needs to reflect the importance and the complexity of these phenomena. In the present model, convective heat losses were simulated by accounting for both free and forced convection and their variation over the human body. Simulation of the radiant heat exchange with the surroundings considered spatial asymmetries of the environment and the (posturedependent) degree of exposure of individual sites of the body elements. Respiration was modeled by considering both the convective and latent heat losses. Shortwave radiation was also included to enable the effects of solar irradiation and other hightemperature sources on the human body to be analyzed. A clothing model that accumulated nonuniform insulation properties was developed for use with the multisegmental model. In the skin evaporation model, the traditional skin wettedness concept was replaced by ensuring local heat and mass balances. These considered moisture diffusion through the skin, the evaporation of sweat, the evaporative resistance of garments, and the accumulation of the liquid sweat on the skin.
A numerically stable finitedifference formulation of the bioheat equation was developed, which applied—with slight modifications—to the steady state and transient heat transport in body tissues that were either cylindrical or spherical. The human body was described by timeindependent material matrices and timedependent blood matrices. Arterial temperatures were solved directly, which was much less computationally demanding than the conventional iteration techniques.
The numerical formulation of the passive system was verified by comparison with known analytic solutions for conduction in cylinders and spheres and for an analytic solution of the bioheat equation in cylinders. These comparisons indicated that the model, with its current geometric discretization, was capable of accurately predicting internal body temperatures and that it was computationally stable, even when unrealistically long time steps were used to model fast transient conditions.
The passive system model thus forms a solid foundation on which an active system can operate. Working together, detailed predictions of the thermal state of human beings, particularly when in manmade enclosures (buildings, vehicles, etc.), should be possible. The complete system and its validation against field measurements, and predictions by the models of others, are the subject of the paper in preparation.
Acknowledgments
The authors express their gratitude to the Deutscher Akademischer Austauschdienst for BritishGerman Academic Research Collaboration; to the British Council; the JosephvonEgle Institut; and to Knödler Stiftung for project funding.
Footnotes

Address for reprint requests and other correspondence: D. Fiala, Institute of Energy and Sustainable Development, De Montfort University, The Gateway, Leicester LE1 9BH, UK.

Submitted 19 July 1996; accepted in final form 24 June 1999.
 Copyright © 1999 the American Physiological Society
Appendix
A finitedifference scheme was used to discretize Eq. 1.The partial derivatives with respect to radius were approximated by using a central difference method, whereas the CrankNicolson approach was adopted for derivatives with time. The CrankNicolson method (Ω_{CN}T_{r}) is constituted by averaging the explicit (Ω_{expl}T_{r}) and the implicit (Ω_{impl}T_{r}) methods, which have been defined for the current time step (t) and the future time step (t + 1), respectively
By applying Eq. EA1
and separating the future temperature terms, one obtains commonly
The time step Δt approximates the differential dt, and β_{r} represents the timedependent energy equivalent of the nodal blood flow rate (see above).
This particular numerical formulation of the bioheat equation has universal applicability. It is valid for heat transfer both in cylindrical and spherical coordinates, and the steadystate expression (used for computing setpoint temperatures in the model) arises by dropping the “transient” coefficient ζ/Δt and by removing all temperature terms on the current (righthand) side
The accuracy of the finite difference method depends on the formulation of the boundary conditions. Usually, temperature nodes are located at interfaces, and the corresponding heat fluxes are calculated assuming a plane geometry of the interfaces. This method, however, may cause errors in transient conditions (32). To avoid this, the last node of a homogenous tissue layer was attached to an additional, imaginary node with the same thermophysical tissue properties. Similarly, an additional, fictitious node was placed before the first node of the next homogenous tissue layer. The boundary condition was then satisfied by equating both the interface temperatures and the heat fluxes between pairs of regular and imaginary nodes of the neighboring material layers (32). Hereby, the impact of the geometry on the heat dissipation in cylinders and spheres was modeled. In cylindrical coordinates, the interface temperature T_{ifc} depends on the interface radiusr
_{ifc} and the temperatures of the adjacent (regular or fictitious) nodes T_{−} (corresponding to the radiusr
_{−} = r
_{ifc} − Δr/2) and T_{+} (corresponding to the radiusr
_{+} = r
_{ifc} + Δr/2)
The boundary condition in the core of each body element was formulated by establishing an isothermal core element around the cylinder axis or the midpoint of the sphere. As an example, Fig. 1 shows this domain for a leg (node no. 1). Because there is no temperature gradient with depth in these isothermal domains, the heat storage element method had to be applied, substituting the conduction term in Eq. 1
by heat fluxes q
_{ifc,j} passing the interface between the core and the adjacent tissue sectors
The conductive heat flux at the surface of a skin sectorq
_{sk}, i.e., the heat flux calculated between the last skin node and the adjacent imaginary node, was calculated fromEqs. EA7
and
EA9
. It is balanced by the sum of the dry heat loss through the clothing and the evaporative heat loss from the skin
The heat transfer coefficient U*_{cl}incorporates three heattransport mechanisms: conduction through clothing, surface convection, and longwave radiation to surrounding surfaces. Because of the way Eq. 24 forU*_{cl} and Eq. 25 for the evaporative coefficient U*_{e,cl} were formulated,Eq. EA11 is valid even when no garments cover the skin. The resultant partial water vapor pressure acting on the skin, i.e., P_{sk} in Eq. EA11, is calculated by Eq. 28.
To obtain the heat dissipation of the whole body, Eq. EA3 was applied to all the tissue nodes of the passive system and was coupled by boundary conditions at the interfaces. The resulting set of linear equations was solved at each time of an exposure for the boundary conditions at the skin surfaces (Eq. EA11 ) and the regulatory responses of the active system (see Ref. 16).
The generated equation set contains both timeindependent coefficients and timedependent coefficients. The timeindependent constants γ_{r}, δ_{r}, and ζ_{r} in Eq. EA3, which describe geometry and thermophysical material properties, have only to be computed once. These terms, taken from the lefthand side of Eq. EA3
, are
The two terms on the lefthand side of Eq. EA3, containing the coefficient
A direct solution for calculating the arterial blood temperatures was developed. This produces numerically exact results and is considerably quicker than the popular but less accurate and computationally expensive iteration techniques. The solution was obtained by substituting elements T_{bl,a} in Eq. EA3
with the numerical equivalent of Eq. 10, in which T_{bl,v} and the integrals had been replaced by Eq. 11, and the corresponding sums, respectively
The blood coefficient matrix of a body element is constituted by collecting the blood flow terms on the lefthand sides of Eq.EA3
but excluding the T_{bl,p} term. Expression A15 shows a simplified blood coefficient matrix [the (t + 1) index was dropped] for a (fictitious) fournode body element
In central body elements where h_{x} is zero (see Table 2), the blood coefficient matrix of a fournode element would reduce to
Equation EA18 completed the set of equations that formed the main matrix of the whole body.
The coefficients of the main coefficient matrix (Fig.8) were generated by adding the timeindependent conduction submatrices and the corresponding timedependent blood submatrices.
The coefficients of the blood pool temperature line in Fig. 8 arise from Eq. EA18
and are
The coefficients of the blood pool temperature columnCc _{j,r} were generated from Eq.EA3 (with Eq. EA14) and result inCc _{j,r} = −Lc _{j,r}δ_{j,r}/ V_{j,r}.
The socalled “coupling coefficient,” which plays a role in the solution procedure, was obtained from Eq. EA18.
A quick solver was developed for the matrix system, which works very efficiently because only the nonzero cells are considered. The solver reduced the computer time by ∼75% when compared with a traditional Gaussian algorithm, and the demand for computer memory was reduced by ∼90%.