Abstract

A dynamic model predicting human thermal responses in cold, cool, neutral, warm, and hot environments is presented in a two-part study. This, the first paper, is concerned with aspects of the passive system:1) modeling the human body, 2) modeling heat-transport 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 well-established 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 heat-transfer 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 heat-transport 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 A-A′ 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 A-A′) 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 4-cm radius.

Subscripts and superscripts

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

Crank-Nicolson

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

Short-wave 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 so-called “bioheat equation.” This differential equation describes the heat dissipation in a homogenous, infinite tissue volumek∂2T∂r2+ωr∂T∂r+qm+ρblwblcbl(Tbl,a−T)=ρc∂T∂tEquation 1

According to Eq. 1, from left to right, the radial heat flow from warmer to colder tissue regions [heat-conduction 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 metabolismqm(W/m3) and blood perfusion [heat-convection term, where ρbl is density of blood (kg/m3), wbl is blood perfusion rate (s−1), cbl is heat capacitance of blood (J ⋅ kg−1 ⋅ K−1), and Tbl,a is arterial blood temperature (°C)]. This combined effect is balanced by the storage of heat within the tissue mass [right-hand side of Eq. 1, where ρ is tissue density (kg/m3), 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 heat-generation term qm, and basal blood perfusion rate wbl for each tissue layer (Table 2). It should be noted that the magnitudes of the physiological variablesqm and wbl are affected by responses of the active system, as described elsewhere (16). The arterial blood temperature Tbl,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 heat-conduction 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 qm andwbl 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 qm,bas,0 and the additional heat Δqm, which may be produced by local autonomic thermoregulation, while the subjects are exercising and/or shiveringqm=qm,bas,0+ΔqmEquation 2

The basal metabolic rates of individual tissue materialsqm,bas,0 were obtained from literature (10, 17, 40) and are listed in Table 2. In thermal neutrality (Tair = 30°C; for a detailed specification see Table3), where no regulation occurs, the overall basal metabolism Mbas,0 = ∫qm,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 Δqm may contain three components; i.e., changes in the basal metabolism (Δqm,bas) and in additional metabolism by shivering and working (qm,sh andqm,w, respectively)Δqm=Δqm,bas+qm,sh+qm,wEquation 3

The change in the basal metabolism Δqm,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 qm,bas,0 is present and where the tissue temperature differs from its setpoint T0. This local regulation arises from the van't Hoff Q10effect with a sensitivity coefficient of 2 (40), which reflects the dependence of biochemical reactions on local tissue temperature TΔqm,bas=qm,bas,0·[2(T−T0)/10−1]Equation 4The shivering term qm,sh in Eq.3 is a portion of the overall regulatory response elicited by the active system (16). The qm,w term can be written asqm,w=∂(am,wH)∂VmusEquation 5where am,w (see Table 2) is a coefficient distributing to body elements the overall extra metabolism due to exercise, Vmus is the corresponding body-element muscle volume, and H is the internal whole body workload. The values of am,w for standing activities were obtained from Stolwijk (34). The present model uses other estimates ofam,w for sedentary activities, because the extra metabolism shifts from the muscles of the legs to the musculature of the trunk and arms when the subject is seated.

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 heat-generation rate remaining in the body and the basal metabolism Mbas,0H=actMbas,0actbas(1−η)−Mbas,0Equation 6

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/m2 (2), i.e., 1 met. In Eq. 6, the whole body metabolism is calculated from the basal metabolic rate Mbas,0 for a reclining, resting individual, corresponding to actbas = 0.8 met. The act units, used as (time-dependent) 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)η=0.2·tanh(b1·act+b0)Equation 7

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 Tbl,a, instead of assuming a constant arterial blood temperature for all body elements equal to the temperature of the central blood pool Tbl,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 Tbl,p − Tbl,a corresponds to the increase in temperature Tbl,v,x − Tbl,v of the blood in the veins of the same body element after CCXρblcbl∫wbldV(Tbl,p−Tbl,a)=ρblcbl∫wbldV(Tbl,v,x−Tbl,v)Equation 8Here, ∫ wbldV represents the blood flow in the veins as the volume-integral of capillary blood flowswbl over the body element. The blood temperatures after undergoing CCX, Tbl,a and Tbl,v,x, are the two unknowns in Eq. 8. To obtain an element's arterial blood temperature Tbl,a in Eq. 8, the net heat exchange between adjacent arteries and veins Qx was defined by using the equation of Gordon (17)Qx=hx(Tbl,a−Tbl,v)Equation 9where hx is the corresponding countercurrent heat exchange coefficient and Tbl,v is the venous blood temperature before CCX. The arterial blood temperature Tbl,a of a body element is then obtained by substituting the right-hand side of Eq. 8 by the right-hand side of Eq.9 and rearrangingTbl,a=ρblcbl∫wbldV·Tbl,p+hx·Tbl,vρblcbl∫wbldV+hxEquation 10

The computation of the central blood pool temperature Tbl,putilized a special calculation process (see numerical methods). Because the bioheat equation assumes that capillary blood reaches an equilibrium with the surrounding tissue, Tbl,v of a body element in Eq. 10 is obtained as followsTbl,v=∫wblTdV∫wbldVEquation 11

The heat-exchange coefficients hx in Eq. 10 are zero for central body elements (17) (i.e., no countercurrent heat exchange). Few reliable data on hx of extremities exist in literature. Therefore, hx values of extremities and shoulders were estimated by a trial-and-error procedure, in which the coefficients were adjusted to obtain agreement between local skin temperatures as predicted and measured. For legs, feet, arms, and hands, hx 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 hx are listed in Table 2.

In thermal neutrality, tissues are supplied at basal perfusion rateswbl,0 (see Table 2). In nonneutral conditions or during work, the blood flows, wbl in Eq. 1,vary with changes in regional metabolic rates qm. Experimental estimations (30) indicate a linear relationship betweenwbl and qm. Instead of dealing with blood perfusion rates wbl, the present model considers directly their energy equivalent β = ρblcblwbl(W ⋅ m−3 ⋅ K−1). The increased demand for oxygen is therefore accounted for by calculating the change in the gain factor Δβ = ρblcblΔwbl as a function of the change in metabolism Δqm using a proportionality constant μbl = 0.932 (K−1), which was obtained from Stolwijk (34)Δβ=μbl·ΔqmEquation 12Utilizing Δβ, the blood perfusion term in Eq. 1becomes β(Tbl,a−T)Equation 13whereβ=β0+ΔβEquation 14

The very variable blood perfusion rates in the cutaneous plexus arise from human temperature-regulation 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 qc, by radiation with surrounding surfacesqR, by irradiation from high-temperature sourcesqsR, and by evaporation of moisture from the skinqe. 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 qsk (W/m2) passing the surface of a peripheral sector is equivalent to the sum of these individual heat exchangesqsk=qc+qR−qsR+qeEquation 15

The elements of this equation are considered in detail in following sections.

Convection.

The convective heat exchange qc between a skin sector of surface temperature Tsf and ambient air of temperature Tair was modeled by considering both natural and forced convection using combined convection coefficients hc,mixqc=hc,mix·(Tsf−Tair)Equation 16

In the present model, the convection coefficients hc,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 airspeedvair,eff (m/s)hc,mix=anatTsf−Tair+afrcvair,eff+amixEquation 17

The basis for hc,mix are the experiments of Wang (36-38), who measured local convective heat losses on a heated full-scale manikin with a realistic skin temperature distribution for both natural and forced convection. Coefficientsanat, afrc, andamix (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 heat-exchange model, the overall, i.e., the mean, convection coefficient of the human body, hc,m(W ⋅ m−2 ⋅ K−1), was obtained by integrating the local convective heat lossesqc = hc,mix(Tsk − Tair) over the (naked) body surface ADuhc,m=∫qcdAskADu(Tsk,m−Tair)Equation 18Here, Tsk,m is the area-weighted mean skin temperature, i.e., Tsk,m = ∫TskdAsk/ADu, where Tsk is the model's local skin temperature, dAsk is the corresponding local surface area at a body element of the humanoid, and ADu = ∫ dAsk is the Dubois area of the passive system, which is 1.86 m2.

The overall coefficient was computed both for conditions of the natural convection (vair,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 2-m-high 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 vair,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 hc,m with increasing body movement. More recent experiments (8), however, have illustrated constant or even decreasing hc 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, hc is affected by the type of activity rather then by its level. Hence, instead of introducing a fixed relationship between hc and the activity level, which could produce systematic errors in predictions, users of the model may modify the input variable vair,eff appropriately, if detailed information on hc for a specific task is available.

Radiation.

The exchange of heat by long-wave radiation qR is usually of similar importance in the heat balance of the human body as the heat exchange by convection.

In asymmetric environments, qR 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 ofqR 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 direction-dependent average temperatures for the surrounding structures was adopted. In the model, either Tsr,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 Tsr,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 Tsr,m, and rearranging the Stefan-Bolzmann law by introducing the (local) radiative heat-exchange coefficient hR(W ⋅ m−2 ⋅ K−1), the energy exchange is calculated byqR=hR·(Tsf−Tsr,m)Equation 19wherehR=ςεsfεsrψsf­sr(T*sf2+T*sr,m2)(T*sf+T*sr,m)Equation 20and ς = 5.67 × 10−8W ⋅ m−2 ⋅ K−4 is the Stefan-Bolzmann constant; εsf and εsr,mare the emission coefficients of the body surface sector considered and of the surrounding surfaces, respectively; ψsf-sr is the corresponding view factor; T*sf and T*sr,m are the absolute temperatures (K) of the body surface sector and of the surrounding surfaces seen by the body sector.

The long-wave 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 Tsr,m, εsr has to be set to unity (black envelope).

The concept of using mean surrounding temperatures Tsr,mmakes it possible to introduce the view factors ψsf-srbetween sectors and the surroundings seen by them. The view factors ψsf-sr of the cylindrical and spherical body elements (Table 2) were calculated—for each sector separately—by a specially developed finite-element 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 body-radiation area. As a check on the reliability of the view factors, the effective radiant area ratio of the (nude) body feff = ∫ ψsf-srdAsk/ADuwas computed. The resultant whole body values offeff = 0.80 for the standing position andfeff = 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 so-called climate-input file. For inhomogeneous environments, Tsr,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 Tsr,m for additional directions. For example, Tsr,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 Tsr,m must be calculated from measurements of individual surface temperatures. A simple approach to calculate Tsr,m is to weight the surface temperatures Tsr, j by the corresponding surface areasAsr, jTsr,m=∑j=1nTsr,jAsr,j∑j=1nAsr,jEquation 21

A more detailed method is to consider the effect of the body's location in the room by using the view factors ψbd-sr, jTsr,m=4∑j=1nψbd­sr,j·T*sr,j4−273.15Equation 22

The overall view factors ψbd-sr,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 high-temperature sources (sun, fireplaces, etc.) was also considered in the formulation of the passive system. The short-wave radiation intensity and its direction is a characteristic of the environment and is provided, together with other environmental variables, in the climate-input file. The amount of heatqsR (W/m2) absorbed at a sector surface is taken into account in the heat balance of superficial body element sectors by the termqsR=αsfψsf­srsEquation 23where αsf is the surface absorption coefficient and depends on the color of the covering material, s (W/m2) is the radiant intensity, and ψsf-sris the view factor between the sector and the surrounding envelope (see Table 2).

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 Icl, the clothing area factor fcl, and the moisture permeability index icl (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 Icl,fcl, and icl for these garments and from the evaporative resistance of the fabric Re,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 heat-transfer coefficientU*cl(W ⋅ m−2 ⋅ K−1) of a j-layered clothing ensemble worn on a body-element sector is computed in the model asU*cl=1∑j=1n(I*cl)j+1f*cl(hc,mix+hR)Equation 24where (I*cl)j(W ⋅ m−2 ⋅ K−1) is the “local” heat resistance of j-th clothing layer,f *cl (dimensionless) is the local clothing area factor of the outer clothing layer, and hc,mix and hR(W ⋅ m−2 ⋅ K−1) are the actual local coefficients for convection and radiation, respectively.

The corresponding evaporative coefficientU*e,cl(W ⋅ m−2 ⋅ Pa−1) is obtained by using the local values i*cland I*cl of individual clothing layers applied to the body sector, f *cl of the outer clothing layer, the local convection coefficient hc,mix, and the Lewis constant for air Lair = 0.0165 K/Pa (24)U*e,cl=Lair∑j=1nI*cli*clj+1f*cl·hc,mixEquation 25

As a check on the validity of the new clothing model, the overall insulation Icl of clothing ensembles and of individual garments was computed by using a similar procedure to that employed for obtaining the overall convection coefficient hc,m. However, ∫qcdAsk,ADu, and Tair in Eq. 18 were replaced by the integral of the local heat fluxes through the clothing ∫ qcldAcl, the overall surface area Acl of the clothed body, and the corresponding mean surface temperature Tcl,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 so-called “KSU-uniform,” Icl= 0.6 clo, icl = 0.34, the overall insulation was predicted to be Icl = 0.59 clo andicl = 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 Ask is given byU*e,cl(Psk−Pair)=λH2OdmswAskdt+Posk,sat−PskRe,skEquation 26

The left-hand side of Eq. 26 represents the net energy transfer as driven by the evaporative potential between skin and air, where Psk is the water vapor pressure at the skin surface, Pair is the vapor pressure of the ambient air, andU*e,cl is the resultant evaporation coefficient of garments covering the sector.

The first right-hand term of Eq. 26 considers the evaporation of sweat from the skin surface, whereby λH2O = 2,256 kJ/kg is the heat of vaporization of water and dmsw/dt is the rate of sweat production over Ask, as elicited by the active system (16).

The last term of Eq. 26 describes the heat transport by moisture diffusion through the skin where Posk,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 Posk,sat depending on the tissue temperature of the outer skin, Tosk (24)Posk,sat=100×exp18.956−4,030Tosk+235Equation 27The model uses a skin moisture permeability equal to 1/Re,sk = 0.003 W ⋅ m−2 ⋅ Pa−1, which is also the published value (14). Simulations using the present model showed that this value produces the model's basal skin wettedness of wtsk = 0.06, which is a basic property of the human skin (2).

The net evaporative heat loss of a skin sector arises from the actual vapor pressure Psk (Pa) found at the skin surface. This quantity is calculated by using a rearranged Eq. 26Psk=λH2OdmswAskdt+Posk,satRe,sk+U*e,clPairU*e,cl+1Re,skEquation 28

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 Psk, as computed from Eq. 28,exceeds the saturated vapor pressure Psk,sat, which is calculated from Tsk by using an equation similar to Eq.27. The rate of moisture storage dmacc/dt(kg/s) is then equivalent to the rate of moisture production dmsw/dt minus the rate of moisture evaporationdmaccdt=dmswdt−U*e,cl(Psk,sat−Pair)λH2OAskEquation 29

According to the experimental results quoted by Jones and Ogawa (20), the maximum amount of sweat liquid on the skin is limited to about macc/Ask = 35 g/m2. 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., ∫qmdV (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 Tair(°C) and the partial vapor pressure of the ambient air, Pair (Pa)Ersp=4.373·∫qmdV(0.028−6.5×10−5Tair−4.91×10−6Pair)Equation 30

Because the enthalpy of the expired air still depends to a certain degree on the condition of the inspired air, Erspappears as a function of Tair and Pair, allowing Ersp 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, Tair and Pair, respectivelyCrsp=1.948×10−3·∫qmdV(32.6−0.066·Tair−1.96×10−4Pair)Equation 31

The total respiratory heat loss Ersp + Crsp 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 arsp 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 heat-generation termqm in Eq. 1qm=qm,bas,0+Δqm−∂[arsp(Ersp+Crsp)]∂VEquation 32

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 finite-difference scheme was used to discretize Eq 1.The partial derivatives with radius were approximated by the central difference method, which provides improved (i.e., second-order) 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 Crank-Nicolson 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 time-independent “conduction” matrices, which collected all thermophysical tissue constants, and time-dependent “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 passive-system 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 warm-air or hot-air 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 steady-state conditions. These tests indicated that the model predicts steady-state 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 qm = β = 0. At t = 0, the sphere, which had an initially homogenous temperature of 37°C, was exposed to ambient temperatures of Tair = 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 = Tair = 0°C. At t = 0, the cylinder is suddenly supplied with a constant rate of blood flow (wbl = 0.001 s−1, Tbl,a = 37°C) and source of metabolic heat (qm = 600 W/m3). 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 time-step length was investigated by using the previous test but with a higher blood perfusion rate (wbl = 0.005 s−1), to ensure an even more rapid heating of the cylinder (reaching ∼95% of the steady-state 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.

Comparison of predicted muscle tissue temperatures in a leg-cylinder with exact, analytically derived values after instantaneous introduction of a high blood flow and metabolic heat. Prediction time step was 10 min.

Finally, the temperature distribution within, and at the surface of, the whole passive system was calculated for conditions of “thermal neutrality” of Tair = 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 Tsk,m = 34.4°C, and a head core temperature (hypothalamus) of Thy = 37.0°C, values that agree with published data of Tsk,m = 34.4°C (14) and Tty = 37.0°C (4) for reclining subjects exposed to thermoneutral conditions. Some further simulation results of this steady-state 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 (posture-dependent) degree of exposure of individual sites of the body elements. Respiration was modeled by considering both the convective and latent heat losses. Short-wave radiation was also included to enable the effects of solar irradiation and other high-temperature 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 finite-difference 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 time-independent material matrices and time-dependent 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 man-made 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 British-German Academic Research Collaboration; to the British Council; the Joseph-von-Egle Institut; and to Knödler Stiftung for project funding.

Appendix

A finite-difference scheme was used to discretize Eq. 1.The partial derivatives with respect to radius were approximated by using a central difference method, whereas the Crank-Nicolson approach was adopted for derivatives with time. The Crank-Nicolson method (ΩCNTr) is constituted by averaging the explicit (ΩexplTr) and the implicit (ΩimplTr) methods, which have been defined for the current time step (t) and the future time step (t + 1), respectivelyΩCNTr(t)=12[ΩexplTr(t)+ΩimplTr(t+1)]Equation A1When applied to the Eq. 1, the explicit and the implicit numerical expressions for the node r of a cylinder yieldΩexplTr(t)=ρrcrTr(t+1)−Tr(t)Δt=krTr+1(t)+Tr−1(t)−2Tr(t)Δr2+Tr+1(t)−Tr−1(t)2rΔr+qm,r(t)+βr(t)[Tbl,a(t)−Tr(t)]ΩimplTr(t+1)=ρrcrTr(t+1)−Tr(t)Δt=krTr+1(t+1)+Tr−1(t+1)−2Tr(t+1)Δr2+Tr+1(t+1)−Tr−1(t+1)2rΔr+qm,r(t+1)+βr(t+1)[Tbl,a(t+1)−Tr(t+1)]Equation A2

The time step Δt approximates the differential dt, and βr represents the time-dependent 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 steady-state expression (used for computing set-point temperatures in the model) arises by dropping the “transient” coefficient ζ/Δt and by removing all temperature terms on the current (right-hand) side(γr−1)Tr−1+(2+δrβr)Tr−(γr+1)Tr+1−δrβrTbl,a=δrqm,rEquation A5

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 Tifc depends on the interface radiusrifc and the temperatures of the adjacent (regular or fictitious) nodes T− (corresponding to the radiusr− = rifc − Δr/2) and T+ (corresponding to the radiusr+ = rifc + Δr/2)Tifc=lnrifcrifc+Δr/2lnrifc−Δr/2rifc+Δr/2T−−lnrifcrifc−Δr/2lnrifc−Δr/2rifc+Δr/2T+Equation A6The corresponding heat flux qifc can be written asqifc=kr(T−−T+)rifclnrifc+Δr/2rifc−Δr/2Equation A7Similarly, Tifc and qifc in a spherical element (head) are given byTifc=rifc+Δr2ΔrT+−rifc−Δr2ΔrT−+T−−T+rifcr−Δr−rifcr+ΔrEquation A8qifc=kr(T−−T+)rifc2rifc−Δr/2−rifc2rifc+Δr/2Equation A9

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 qifc,j passing the interface between the core and the adjacent tissue sectorsρ1c1dT1dt=qm,1+β1(Tbl,a−T1)−κπr1∑j=1sectorsϕjqifc,jEquation A10In Eq. EA10, the symbol ϕj is the angle of sector j (see Fig. 1 for ϕ j of the anterior leg), and it arises—together with the factor κ/πr1—from relating the interface area of the jth sector to the core volume. The parameter κ depends on the geometry and is κ = 1 for cylinders and κ = 3/2 for spheres. The conductive heat fluxesqifc,j are obtained by Eqs. EA7 and EA9 for cylinders and spheres, respectively. The Crank-Nicolson scheme was also used to discretize Eq. EA10.

The conductive heat flux at the surface of a skin sectorqsk, 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 skinqsk=U*cl(Tsk−To)+U*e,cl(Psk−Pair)Equation A11where the operative environmental temperature To(°C) was formulated to integrate the influences of convection (hc,mixTair), the temperature radiation (hRTsr,m), and short-wave radiation (ψsf-srαsfs) absorbed by a body sector surfaceTo=hc,mixTair+hRTsr,m+ψsf­srαsfshc,mix+hREquation A12

The heat transfer coefficient U*clincorporates three heat-transport mechanisms: conduction through clothing, surface convection, and long-wave 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., Psk 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 time-independent coefficients and time-dependent coefficients. The time-independent 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 left-hand side of Eq. EA3, are(γr−1)Tr−1(t+1)+ζrΔt+2Tr(t+1)−(1+γr)Tr+1(t+1)Equation A13They were arranged to form the tridiagonal conduction matrices, which are a characteristic of the finite-difference method. The pattern for the leg of the body (see Fig. 1) is shown in Fig.7, with the off-diagonal terms reflecting the coupling role of the core element.

The two terms on the left-hand side of Eq. EA3, containing the coefficient βr(t+1), formed “blood matrices” for each body element. Because the coefficients of a blood matrix change with time, they had to be determined for each time step and each step of any iteration procedure simulating the active system.

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 Tbl,a in Eq. EA3 with the numerical equivalent of Eq. 10, in which Tbl,v and the integrals had been replaced by Eq. 11, and the corresponding sums, respectivelyTbl,a=Tbl,p∑rnodesβrVrhx+∑rnodesβrVr+hx∑rnodesTrβrVr∑rnodesβrVrhx+∑rnodesβrVrEquation A14

The blood coefficient matrix of a body element is constituted by collecting the blood flow terms on the left-hand sides of Eq.EA3 but excluding the Tbl,p term. Expression A15 shows a simplified blood coefficient matrix [the (t + 1) index was dropped] for a (fictitious) four-node body elementδ1β1+Bxβ12Bxβ1β2Bxβ1β3Bxβ1β4Bxβ2β1δ2β2+Bxβ22Bxβ2β3Bxβ2β4Bxβ3β1Bxβ3β2δ3β3+Bxβ32Bxβ3β4Bxβ4β1Bxβ4β2Bxβ4β3δ4β4+Bxβ42Equation A15whereBx=−hxδrVr∑rnodesβrVrhx+∑rnodesβrVrEquation A16

In central body elements where hx is zero (see Table 2), the blood coefficient matrix of a four-node element would reduce toδ1β10000δ2β20000δ3β30000δ4β4Equation A17Because the central blood pool temperature Tbl,pappears (as a further unknown) in Eq. EA14 for all the nodal heat balances, an extra equation was needed that balances the heat content of the pool's inlet and outlet. The equilibrium blood pool temperature resulted in an expression similar to Eq. 11;however, the product of the elements' venous blood temperatures after CCX, Tbl,v,x, and the corresponding blood flow rates, as well as the local perfusion rates wbl, were integrated over the whole body. In the final Eq. EA18,Tbl,p is a function of all the nodal tissue temperatures and underlines the coupling effect of blood circulation in the heat dissipation of the passive systemTbl,p=∑jelem∑rnodesβj,rVj,rhx,j+∑rnodesβj,rVj,r×∑rnodesβj,rVj,rTj,r∑jelem∑rnodesβj,rVj,r2hx,j+∑rnodesβj,rVj,rEquation A18

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 time-independent conduction submatrices and the corresponding time-dependent blood submatrices.

The coefficients of the blood pool temperature line in Fig. 8 arise from Eq. EA18 and areLcj,r=∑rnodesβj,rVj,rhx,j+∑rnodesβj,rVj,r×βj,rVj,rEquation A19

The coefficients of the blood pool temperature columnCcj,r were generated from Eq.EA3 (with Eq. EA14) and result inCcj,r = −Lcj,rδj,r/ Vj,r.

The so-called “coupling coefficient,” which plays a role in the solution procedure, was obtained from Eq. EA18.CplC=∑jelem∑rnodesβj,r2hx,j+∑rnodesβj,rEquation A20

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%.